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I.  Introduction 

Problems  of  hyperveloclty  Impact  attract  attention  of 
scientists  during  the  last  thirty  years  and  occupy  a  significant 
part  In  scientific  literature.  The  urgency  of  Investigations  of 
these  problems  follows  from  the  practice.  It  concerns,  for 
example,  evaluation  of  effectiveness  of  spacecraft  debris 
shields,  analysis  of  asteroid  Impact  events,  problems  of  atomic 
power  stations  safety  and  so  on.  From  the  other  side  the 
possibility  of  numerical  simulations  of  these  phenomena  becomes 
real  owing  to  progress  In  theoretical  and  experimental  high 
pressure  high  temperature  physics.  The  success  of  numerical 
simulation  of  hypervelocity  Impact  depends  on  two  factors: 
physical  model  (Including  mechanical,  kinetic, thermodynamic, 
etc.  properties  of  matter)  and  numerical  methods  for  Integration 
of  governing  equations.  The  second,  part  of  the  problem  will  be 
discussed  below.  As  to  the  flrst^one,  the  review  of  various 
substances  models  Is  given,  for  example.  In  [1,  23. 

The  most  widespread  numerical  methods  for  tiyperveloclty 
Impact  (HVI)  simulation  are  based  on  Langranglan  and  particle 
techniques  [  3-7  3 .  Huge  experience  has  been  stored  up  to  now  In 
this  area.  It  becomes  clear,  that  the  main  drawback  of 
Lagranglan  codes  Is  the  restriction  on  deformations  of'  flow  and 
that  of  for  the  PIC  codes  Is  the  necessity  of  keeping  of 
’’superfluous”  Information  about  flow  parameters  In  additional 
numerical  grid  cells  Into  which  substance  can  get.  Of  course, 
these  difficulties  are  algorlt3imlc  and  can  be  overcome  by  some 
additional  measures.  For  example,  a  usual  procedure  for 
Lagranglan  teclmlques  Is  re  Interpolation  of  flow  parameters 
into  a  new  grid,  when  the  grid  deformations  become  large  183.  As 
to  Eulerlan  tec3inlques,  one  can  mention  about  tracking  methods- 
[93,  that  are  used  to  Increase  the  accuracy  of  calculations  near 
matter-vacuum  and  contact  boundaries  and  to  restrict  the  domain 
of  computations.  Unfortunately  Eulerlan  codes  have  .one  more 
disadvantage  associated  with  mass  dispersion  problems,  that 
makes  It  difficult  to  calculate  debris  cloud  evolution.  The  are 
several  ways  to  decrease  this  dispersion.  For  example  In  [103 
both  the  projectile  and  the  shield  were  proposed  to  have  nonzero 
Initial  velocities  moving  towards  each  other  to  minimize 


advectlon  ol  materials  between  cells  in  the  Eulerlan  mesh.  In 
»  order  to  decrease  dispersion  In  Eulerlan  computations  hlgh-order 
accurate  advectlon  schemes  also  were  used  [11,  12]. 

As  both  pure  Lagranglan  and  Eulerlan  codes  have  their  own 
drawbacks,  It  Is  naturally  to  develop  arbitrary  Lagrange-Euler 
(ALE)  codes  combining  the  advantages  of  these  approaches. 
Following  this  Idea,  numerical  simulation  of  penetration  of  a 
hypervelocity  projectile  through  multi-plate  shield  was  done 
[13]  using  two  codes:  one  for  computing  of  penetration  stage 
(Eulerlan)  and  another  (Lagranglan)  for  debris  cloud  transport 
calculation.  In  recent  years,  some  algorithms  of  computation 
of  real  matter  dynamics  (Including  elasto-plastlc  and  damaged 
matter  )  developed  previously  only  for  Lagranglan  frame  have 
been  extended  to  the  Eulerlan  frame  of  reference  [12,  14,  15]. 
This  opens  new  possibilities  for  HVI  simulations  using  more 
flexible  ALE  codes. 

In  the  present  work  we  used  Godunov’s  scheme  for 
computations  [16].  Godunov’s  me  till'd  In  moving  grids  provides 
great  possibilities  and  Includes  both  Lagranglan  and  Eulerlan 
approaches  as  partial  cases.  It  should  be  marked  that  this 
scheme  has  been  Intensively  explored  for  simulations  of 
aerodynamics  problems ,  where  high  accuracy  of  calculations  Is 
necessary. 

In  1968  S.K. Godunov,  A. V. Zabrodin  and  G.P.Prokopov  [IT] 
formulated  a  programmatic  approach  to  numerical  Integration  of 
nonsteady  Euler  equations  on  moving  grids.  During  several  years 
the  authors  of  this  report  have  being  working  on  realization  of 
this  approach.  Results  of  this  activity  are  resented  below. 

For  simplicity  the  consideration  Is  given  for  the  Euler 
equations  omitting  of  strength  effects.  Applications  of  this 
approach  to  the  problems  of  high  velocity  Impact  Illustrate 
possibilities  of  2D  code  developed.  As  distinct  from  the  most- 
Illustrations  of  possibilities  of  this  method  given  in  [17],  the 
problems  of  hyperveioclty  Impact  are  characterized  by  strong 
deformations  of  matter  and  wide  range  of  matter  state 
parameters:  from  Initially  shock-compressed  or  strong^ly  heated 
condensed  state  to  finally  expanding  gaseous  or  spalled 
substance.  Large  variations  of  density  lead  to  an  enormous 
Increase  of  geometrical  sizes  of  computed  fields  occupied  by 
the  substance.  Contacts  separating  the  substance  from  vacuum  are 


the  physical  boundaries  of  computational  region.  Therefore  the 
choice  of  Lagranglan  meshes  seems  to  be  a  preferable  strategy 
for  simulations  of  this  class  of  problems.  Unfortunately 
Lagranglan  meshes  crash,  when  the  deformations  become  large, 
because  the  cells  deform,  twist  and  finally  get  tangled.  The 
first  factor  leads  to  great  restrictions  on  calculation  time 
step  and  the  second  makes  further  calculations  Impossible.  As  to 
the  Eulerlan  grid  the  marked  above  uncertainty  In  computational 
region  leads  to  the  necessity  of  keeping  additional  grid  sells 
and  to  corresponding  Increase  of  computational  time.  The  other 
problem  of  Eulerlan  techniques  Is  the  problem  of  contacts 
between  substances. 

Methods  based  on  the  use  of  computational  regions  with 
moving  boundaries  and  generation  Inside  of  these  regions  of 
grids  connected  only  with  the  boundaries  seems  to  unify 
the  advantages  of  Eulerlan  and  Lagranglan  meshes  allowing  to 
overcome  their  drawbacks.  In  the  most  examples  given  below  the 
emphasis  Is  placed  on  demonstration  of  possibilities  of 
application  of  moving  grids  for  treating  the  flows  with  strong 
deformations.  Thus,  for  example,  simulations  of  Impact  problems 
presented  In  this  report  have  been  performed  on  moving  grids 
with  fitting  all  contacts  between  different  media  Including  the 
boundaries  with  vacuum.  During  the  process  of  Interaction  the 
structure  and  configuration  of  boundaries  change:  segments  of 
contacts  can  disappear, for  example  because  of  rebound,  and  new 
segments  can  originate  due  to  Involving  of  new  bodies  Into  the 
process  of  Interaction.  The  flow  region  under  consideration  can 
be  subdivided  Into  several  subregions.  Boundaries  of  subregions 
can  be  both  as  original  flow  singularities  (characteristics, 
boundary  between  substance  and  vacuum,  the  contact  boundary 
between  two  substances,  the  axis  (plane)  of  symmetry,  the  shock 
wave)  and  so  moving  In  accordance  with  their  physical  nature.- 
They  can  be  also  Inserted  for  some  computational  reasons. 

As  It  was  mentioned  above,  the  basic  attention  Is  paid  to 
computational  aspects  of  simulation.  Nevertheless,  the  algorithm 
developed  makes  It  possible  using  of  different  physical  models 
In  the  frame  of  Godunov’s  scheme.  It  concerns,  for  example,  the 
Rlemann  problem  solver,  designed  for  the  case  of  arbitrary  EOS. 


II.  Physical  model 


Governing  Equations 

One  fluid  one  temperature  hydrodynamic  model  Is  used  In 
simulations.  Governing  Euler  equations  are  written  In  the  form 
of  conservation  laws  of  mass,  momentum,  energy  and  porous 
volume,  what  Is  necessary  for  Godunov’s  method.  The  conservative 
form  of  the  motion  equations  allows  to  describe  both  as 
continuous  and  so  discontinues  flows. 

^  +  vpu  =0, 

+vpu  X  u=  -vp(1-Vp), 

+  vpuh=  -vpud-Vp)^  (2.1) 

+  vpuV^=  fxf. 

Here  p-  density,  u-  mass  velop<ity,  h-total  specific  energy, 
p-  pressure,  V^-  porous  volume  per  unit  volume,  (p-  porous 
growth  rate.  For  p  we  use  semlemplrlcal  spall  damage  kinetic 
built  on  the  base  of  free  surface  motion  registration  [1]. 
Nevertheless  this  model  of  spall  fracture  Is  valid  only  for 
small  (up  to  several  percents)  porosity.  When  the  porosity  of 
matter  becomes  greater  then  some  value,  the  substance  does  not 
resist  to  expansion.  In  this  case  the  volume  occupied  by  the 
solid  component  are  not  changed  anymore  and  the  porosity  grows 
with  the  Increase  of  specific  volume  of  matter. 

Equations  of  state 

Semlemplrlcal  wide-range  equations  of  state  (EOS)  are  used 
[18]  to  determine  the  pressure,  the  temperature  and  the  sound  ' 
velocity  as  functions  of  Internal  energy  per  unit  mass  s  and 
the  density  of  solid  component  p^  connected  with  density  p  as 
Pp=p/(I-Vp).  These  equations  are  built  on  the  base  of 
available  Hugonlot  data  and  provide  for  correct  asymptotes  In 
the  cases  of  extremely  high  energies  and  small  densities. 


III.  Algorithm  of  computations 

1 .  Scheme  of  integration  on  the  moving  grid 

Eqs.(2.1)  are  Integrated  on  the  moving  grid  by  Godunov’s 
procedure  [191.  The  scheme  of  the  algorithm  Is  shown  In  Fig.  1. 
In  comparison  with  the  traditional  one  ([17])  this  algorithm 
enables  a  possibility  of  the  numerical  region  division  Into  many 
subregions  with  Independent  grid  generation  In  every  subregion. 
An  Important  problem  which  should  be  solved  by  that  Is 
Interaction  between  the  boundaries  belonged  to  different 
subregions. 

2.  Movement  of  boundaries 

The  first  step  of  the  algorithm  Is  the  displacement  of 
subregion’s  boundaries.  The  boundary  of  each  subregion  Is 
represented  as  a  coordinate  arr^^  of  nodes.  Any  boundary 
between  adjacent  subregions  represents  a  unification  of  nodes 
belonging  to  boundaries  of  these  subregions.  The  coordinates  of 
nodes  are  also  stored  In  special  arrays  establishing  the 
sequence  of  nodes.  Nodes  corresponding  to  the  beginning  and  to 
the  end  of  different  segments  of  boundaries  are  also  marked.  It 
means  that  each  boundary  segment  Is  represented  as  a  sequence  of 
points  belonging  to  the  boundaries  of  different  subregions  and 
that  of  special  end-points  of  boundaries.  The  total  number  of 
nodes  representing  boundary  exceeds,  generally  speaking,  the 
number  of  cells  being  In  contact  from  the  side  of  adjacent 
subregions. 

To  determine  velocity  of  a  boundary  rib,  the  Rlemann 
problem  with  flow  parameters  from  adjacent  cells  being  In 
contact  along  every  link  Is  solved.  In  case  the  Rlemann  solver- 
yields  a  negative  pressure  a  rebound  along  the  link  Is  supposed 
to  be  happened.  This  fact  Is  taken  Into  account  by  Inserting 
some  new  special  boundary  points  In  the  array.  In  more  complex 
physical  model  the  coupling  force  can  be  prescribed  to  the 
contact  and  a  rebound  can  be  Inserted  when  the  tense  strengths 
exceeds  the  value  of  this  force. 


Flg.1  Block-scheme  of  the  algorithm 


The  velocities  of  nodes  can  he  calculated  after  solving  of 
Rlemann  problem  at  the  rib  connecting  the  nodes.  In  the 
Godunov’s  method  the  flow  parameters;  velocity,  density, energy 
are. assumed  to  be  uniform  Inside  the  meshes.  For  example,  one 
should  take  flow  parameters  from  the  meshes  A  and'  B  to  calculate 
velocity  and  so  from  the  meshes  A  and  C  to  calculate  velocity 
U^,  see  fig. 2  There  may  be  different  ways  of  nodes  velocities 
Interpolation.  The  linear  Interpolation  proposed  by  Godunov  ; 
[17]  gives  the  formula 


12+21 

u=  - — -  (  3.1) 


How  It  ’  s  seen  from  (  3.1)  the  shorter  the  rib  Is ,  the  more 
contribution  gives  It  In  the  node  velocity.  (For  example.  If 
ij»i2  then  u  =  u^)  From  another  .^polnt  of  view.  If  we  want  to 
minimize  the  middle-square  deviation  between  the  boundary 
displacement  obtained  from  the  Rlemann  problem  solution,  when 
every  rib  moves  with  an  appropriate  ”blg”  velocity,  and  the 
boundary  displacement  determined  by  the  nodes  movement  one 
should  use  the  following  formula  (how  It  has  been  shown  In 
[201) 

u  1  1 

*11+22 

u  =  -  (3.2)  , 

1+1 

12 

In  present  algorithm  We  use  an  average  way  of  Interpolation 
between  (3.1  )  and  (3.2  )  ,  which  leads  to  a  simple  formula 

u=0.5*(u  +u*  )=0.5>k  (u^+u^  )  (3.3  ) 

The  formula  (3.3)  combines  the  advantage  of  (3.3)  with 
Simplicity.  ^  ' 


REGION  2 


■.h 

Fig. 2  Boundary  between  the  dl'flerent  subregions  (solid  line) 
’  ^2  velocities  obtained  Irom  the  Rlemann 

problem  solution  between  the  meshes  A  (region  1)  and  B 
(region  2)  and  so  A  and  C,  1^,  1^  -  the  lengths  of 
the  ribs 


3.  Boundary's  type  checking 

Using  of  moving  grid  makes  It  much  more  easier  setting  of 
boundary  condition  In  comparison  with  Elerlan  grid.  In 
accordance  with  the  boundary  condition  some  type  of  boundary  has 
been  assigned  to  every  rib  belonged  to  the  boundary.  There  are- 
five  boundary  types  that  have  been  considered:  "free  boundary” 
(boundary  with  vacuum),  ’’contact”  (boundary  with  another 
matter),  ’’axes”  (axes  of  symmetry),  ’’hard  wall”  and  so  called 
’’special  boundary”  through  which  a  matter  can  flow.  I1^’s  clear, 
that  the  type  of  the  boundary  can  be  changed  from  one  time  step 
to  another.  Therefore  a  procedure  checking  the  boundary  type  at 
each  time  step  has  been  developed. 


4.  Nodes  distribution  along  the  boundary  and  grid 
generation  technique 

Experience  of  the  numerical  simulation  on  the  moving  grid 
shows  that  some  procedure  of  the  nodes  distribution  along  the 
boundary  Is  required  to  generate  a  good  grid  Inside  of  the 
region.  Interpolation  of  the  boundary  plays  an  Important  role 
for  displacement  of  nodes  along  the  boundary.  How  It  has  been 
noticed  [19],  the  sharp  angels  of  the  boundary  are  smoothed 
from  one  time  step  to  another,  when  the  nodes  move  along  the 
stral^t  lines  connecting  them.  To  prevent  this  smoothing  on  the 
boundary,  we  approximate  It  by  arcs  linking  the  boundary  nodes. 

In  the  case  of  strong  deformations  of  numerical  region  the 
law  of  nodes  distribution  along  the  boundary  Influences  In  a 
great  power  on  the  grid  generated  Inside  the  region.  Previously 
[20]  the  law  of  the  nodes  distribution  along  the  boundary  was 
chosen  anytime  In  accordance  with  the  specifics  of  the  problem 
to  be  solved.  Such  a  way  required  ^t  least  several  attempts  to 
find  the  appropriate  law  of  nodes  distribution  for  problem 
considered. 


Fig. 3  The  scheme  ol  the  grid  mapping 


In  this  version  of  moving  grid  code  a  procedure  of  boundary 
node  location  correction  Is  proposed,  which  allows  to  generate 
orthogonal  grids  In  subregions.  For  this  purpose,  boundary  nodes 
are  displaced  along  the  boundary  at  any  time  step  to  ensure 
orthogonality  between  the  boundary  and  correspondent  Inner  grid 
line.  Nodes  shift  Is  limited  by  settled  maximum  ratio  between 
adjacent  the  node  rib’s  lengths.  In  our  calculations  the 
maximum  ratio  was  1.1.  In  order  to  not  distort  the  shape  of 
numerical  region  during  nodes  distribution  the  node  Is  shifted 
along  smooth  curve  line  which  approximates  the  boundary.  The 
curvature  of  the  line  Is  changed  from  node  to  node  using 
line  Interpolation  formula  for  curvature  radius.  After  nodes 
displacement  new  grid  Is  generated  In  subregions. 

Grid  generation  algorithm  Is  based  on  a  univalent  mapping 
of  a  rectangular  region  at  a  parametric  plane  onto  a 
subregion  under  consideration  at  physical  plane  (x,y)  of  flow. 
The  problem  can  be  formulated  mathematically  In  a  following  way. 
The  flow  region  at  plane  (x,y)  Is  .^nsldered  to  be  a  curvilinear 
quadrangle  with  tops  A,  B,  C,  D  (flg.3).  The  univalent  mapping 
of  this  quadrangle  onto  the  rectangular  A',B',C',D'  Is 
determined  by  solution  of  a  system  of  elliptical  equations 

-  ^xx'^^yy^*^  ”^xx'''”'^yy^^  (3.4) 


It  can  be  easily  proved  that  this  mapping  Is  univalent, 
for  example  considering  the  principle  of  maximum  of  harmonic 
function.  In  reality,  to  solve  the  system  of  equations 
relative  to  x,  y  as  Independent  variable  parameters  Is  more 
convenient 

a  x^^  -  2b  x^^  +  g  X^^  =0  (3.5) 

ay^^  -  2b  y^^.  g  y^  =0 
a  =  x^  +  y^  b  =  x^x^  +  y^y^  g  =  x^  +  y^ 


using  the 


following  boundary  conditions: 


x(e)=x(e,'r);L) 

X(T))=X(e(j,Tl) 

x(e)=x(e,T)p) 

x(r|)=x(^^,T]) 


y(U=y(e,T)i) 

y(Ti)=y(e(i,T)) 

y(e)=y(6,Pp) 

y(T))=y(eu,'n) 


at  AB  (3.§) 
at  BC 
at  CD 
at  DA 


where 


Til  ^  “n  <  \ 


The  transition  to  discrete  analog  of  boundary  conditions  for 
this  system  of  quaslllnear  elliptical  equations  can  be 
performed  easily.  If  the  the  grid  nodes  are  given 


1^=1  1=1,..., n  (3.7) 

•  ti~i » •  •  • 

then  the  boundary  conditions  are  (3.8) 

x(l)=x(l,1)  y(l)=y(l,1)  at  AB 

x(J)=x(1,d)  y(l)=y(1,J)  at  BC 

x(l)=x(l,m)  y(l)=y(l,m)  at  CD 

x(J)=x(n,d)  y(d)=y(n,:I)  at  DA 


Different  algorithms  have  been  proposed  In  [IT]  for 
solving  this  problem.  Iterations  considered  to  be  the  most 
efficient  algorithm.  We  modified  this  algorithm  treating 
solution  of  corresponding  equatlons_  as  a  result  of  relaxation 
of  solutions  of  nonsteady  equations.  For  this  purpose  ’’non¬ 
steady”  terms  dx/dt  and  dy/dt  are 'inserted  Into  the  right  sides 
of  equations  for  x  and  y  correspondingly  [21].  The  resulted 
parabolic  system  Is  solved  by  the  method  of  variable  directions 
At  the  first  half  time  step  terms  containing  derivatives  over  h 
are  considered  to  be  given.  At  the  second  half  time  step 
derivatives  over  x  are  taken  from  calculations  at  the  first 
stage.  The  problem  at  the  each  stage  Is  similar  to  the 
one -dimensional  heat  transport  problem.  This  consideration  Is 
similar  to  Iterations,  but  as  any  relaxation  method  It  provides 
more  possibilities  for  the  choice  of  ’’time-step”  or,  what  Is  the 
same,  of  relaxation  parameter.  In  particular.  It  can  be  changed 
from  one  grid  node  to  another  to  accelerate  the  convergence  of 
the  method.  The  proposed  algorithm  has  demonstrated  Its 
reliability.  Nevertheless  It  falls  for  boundaries  of  a- 
complex  configurations.  In  order  to  overcome  problems  associated 
with  a  grid  generation  the  algorithm  foresees  the  possibility  of 
subdividing  the  flow  region  Into  subregions  with  reasonably 
regular  boundaries.  ^ 


5.  Riermnn  problem,  solution  and  flow  parameters 
recalculation  for  next  time  level 

The  basic  point  of  numerical  Integration  of  governing 
equations  Is  the  Rlemann  problem,  which  Is  not  self  similar  If 
any  kinetics  phenomena  Is  considered,  (for  example,  stress 
relaxation  or  mlcroporous  growth).  In  the  model  used  the  porous 
growth  Is  also  Included  Into  consideration.  Therefore  we  take  It 
to  be  "frozen"  when  the  time  of  the  fracture  kinetics  t^^  Is 
much  more  than  the  time  step  At  and  on  the  contrary.  Hugonlot 
adlabates  and  Isoentrops  for  solid  matter  are  taken  for  solving 
Rlemann  problem  In  the  first  case  and  those  for  porous  matter 
are  used  In  the  second  case.  Thls^-approach  has  been  previously 
used  [22]  In  calculations  of  vlsco-elastlc  flows  by  Godunovs 
method. 
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Fig. 4  Discontinuity  decay:  1  and  2  designate  undisturbed 
regions,  1 ’  and  2/  -  regions  behind  the  wayes, 
dashed  line  -  contact  dlscont-lnulty 
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Fig. 5  U-P  plane.  Points  1  and  2  designate  Initial  How 
parameters  from  the  lef.t,i''and  from  the  right  side 
of  the  rib.  U  means  the' normal  to  rib  velocity. 

Unlike  the  perfect  gas,  In  the  case  of  an  arbitrary  EOS 
the  Rlemann  problem  can’t  be  solved  exactly  by  Iteration  scheme 
without  numerous  EOS  evaluations.  One  should  mention  here  about 
some  approximate  Rlemann  solvers  [23,24]  that  can  be  useful  also 
for  the  problems  considered.  Initial  flow  discontinuity  decays 
Into  the  contact  discontinuity,  the  left  wave  and  the  right 
wave,  those  may  be  the  shock  wave  or  rarefaction  wave  depending 
on  the  Initial  parameters  (see  Flg.4).  The  mass  velocities  and 
the  pressures  behind  the  right  wave  (2-2 ’ )  are  equal  those 
behind  the  left  wave  (1-1’).  These  so  called  ”blg  parameters” 
(PG  and  UG)  can  be  determined  as  Intersection  of  left  and  right 
wave  velocities  functions  In  U-P  plane  (Fig. 5  ). 

Densities  and  energies  In  the  regions  1’  and  2’  can  be 
calculated  If  velocity  UG  and  pressure  PG  are  known.  Velocity 
functions  FI  and  FR  are  built  using  a  real  equation  of  state.  To 
build  the  simple  wave  velocity  function  (2 ’-2)  sound  velocity  as 
a  function  of  pressure  and  density  Is  sufficient.  In  the  case  of 
the  shock  wave  (1-1’)  the  energy  as  a  function  of  pressure  and 
density  Is  necessary.  In  order  to  calculate  the  fluxes  of  mass, 
energy,  porous  volume  and  momentum  the  rib  velocity  has  to  be 


known.  Then,  the  fluxes  can  be  determined  choosing  the  flow 
parameter  from  the  region  where  the  rib  velocity  Is  located 
(region  2’  In  the  case  shown  In  Fig. 4  ).  The  velocity  of  the  rib 
V  can  be  determined  If  the  locations  of  the  rib  at  the 
neighboring  time  levels  are  known  (see  Fig. 6).  Generally 
speaking,  the  surface  covered  by  the  rib  during  It’s  movement  Is 
not  the  plane.  For  simplicity  we  replace  It  by  the  plane 
following  the  way  proposed  by  Godunov  [17].  In  accordance  with 
this  method  the  normal  component  of  the  rib  velocity  Is 
calculated  aS 

.^abb’a’ 

V= - 7 — ^  (3.9) 

T  t.  +  A  t  /  2 

*^ab 

Here  A  -  Is  the  area  covered  by  the  rib  ab  during  the  time 
step  At  In  the  plane  XY,  1^^,-the  length  of  the  rib  at  the  time 
t+At/2. 
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Fig. 6  Determination  of  rib  velocity 


6.  Conclusion  remark 

Naturally  the  description  of  the  algorithm  given  above  Is 
at  least  schematic  and  except  of  some  modifications  and  details 
represents  only  the  Interpretation  of  Ideas  [17].  However  we 
considered  the  program  approach  as  being  already  given  and  the 
main  problem  Is  to  put  It  Into  practice  for  specific  problems  of 
physics.  On  this  way  one  can  lace  with  a  lot  of  concrete 
problems  of  various  extent  of  complexity:  mathematical, 
physical,  algorithmic.  One  of  the  most  Important  points  Is  an 
experience  of  working  with  given  complex  of  program.  We  are  sure 
that  the  difficulties  connected  with  realization  of  the  approach 
will  be  rewarded  by  Its  huge  possibilities.  The  method  allows  at 
least  In  principle  to  fit  flow  singularities  :  shocks  , 
contacts,  characteristics.  From  this  point  of  view  the  method 
unifies  main  advantages  of  characteristics  methods  with  the 
advantages  of  uniform  computational  schemes.  A  promising  future 
we  see  In  Integration  of  this  metAod  with  analytical  methods  of 
gas  dynamics.  This  method  provides  a  possibility  to  determine 
dependence  and  Influence  domains  and  to  perform  calculations  at 
this  minimal  domains  area.  From  this  point  of  view  It  can  be 
called  economic. 

Last  years  are  being  marked  by  the  Intense  activity  In 
developing  of  high  order  accuracy  modifications  of  Godunov’s 
methods.  Review  of  these  methods  can  be  found  for  example  In 
[251.  These  modifications  as  a  rule  can  easily  be  Included  Into 
•the  body  of  the  original  Godunov ’ s  algorithm.  Modifications 
concern  the  Rlemann  problem  procedure  and  formulae  for 
recalculation  of  flow  parameters  for  the  new  time  step. 
Development  of  high  order  resolution  schemes  In  moving  grids 
fitting  flow  singularities  Is  of  great  Importance  [261.  They  will 
be  considered  more  detally  In  the  next  report. 


IV.  Demonstrative  examples. 


1 .  TEST  )f9  1 :  impact  of  a  plane  alminm  plate  onto  a 
semiinfinite  aluminum  target 

The  simplest  test  problem  for  2D  code  Is  an  Impact  of  a 
plane  projectile  on  a  semllnflnlte  target.  In  this  case  the 
problem  becomes  to  be  one -dimensional  and  If  Hugonlot  adlabates 
of  materials  are  Known  the  pressure  and  the  shock  velocity  can 
be  easily  calculated.  The  time  when  rarefaction  wave  from  the 
back  side  of  the  projectile  reaches  the  shock  front  can  be  also 
determined.  Results  of  a  plane  aluminum-  on-  aluminum  Impact 
with  8  km/s  velocity  are  shown  In  Fig.  7-8.  Region  of 
computations  consists  of  two  subregions  corresponding  to  the 
target  and  the  Impactor.  Left  boundary  of  the  target  moves 
Initially  with  a  characteristic  velocity  (sound  speed)  Inside 
aluminum  target.  Then  shock  from  contact  reaches  this  boundary 
and  It  becomes  moving  as  shock  dlsj^bntlnulty  (see  fig.  8).  This 
transition  from  sound  to  shock  regime  of  boundary  movement  Is 
automatic  because  the  velocity  of  the  boundary  Is  obtain  from 
Rlemann  problem  solution  algorithm.  As  It  can  be  easily  proved, 
numerical  result  Is  agreed  well  with  the  analytical  calculation 
of  the  shock  pressure  and  velocity. 


2.  TEST  2:  influence  of  the  grid  refinement  on  the 
results  of  numerical  simulation  of  penetration 
problems 

Calculation  with  different  levels  of  accuracy  were 
performed  for  aluminum  sphere  of  1  cm  diameter  penetrating 
aluminum  plate  of  2  mm  thickness  at  10.1  km/s.  Numerical  grids, 
and  levels  of  constant  pressure  corresponding  the  time  moment  of 
710  ns  are  shown  In  flg.9-11.  In  the  case  presented  In  fig.  11 
accuracy  was  two  times  more  than  that  In  the  fig.  10  and  four 
times  more  than  In  the  fig. 9.  It  Is  seen  that  the  higher  the 
accuracy,  the  better  Is  shock  front  resolution  In  the 
projectile.  More  fine  grid  allows  also  to  resolve  the  low 
density  vaporized  matter  blowing  off  the  target  during  impact. 
At  the  following  stage  of  penetration  shown  In  fig.  12-1 3  fine 


grid  gives  less  tMckness  of  the  target  ’  s  substance 

surrounded  debris  of  the  projectile.  And  of  course  as  a 
consequence  of  higher  pressure  there  Is  a  higher  expansion 
velocity  of  debris  cloud  In  this  case. 

3.  Interaction  of  spherical  projectile  with  two  layered 
alminm  target 

Result  of  computation  of  spherical  aluminum  projectile 
Interaction  with  a  spaced  two  layered  aluminum  target  Is 
presented  In  fig. 14.  Impact  velocity  for  the  problem  considered 
Is  8  km/s.  One  should  note  that  not  more  than  1000  meshes  Is 
used  In  this  computation.  The  number  of  meshes  grows  during 
computation  to  ensure  some  fixed  level  of  accuracy.  Criteria  for 
grid  refinement  Is  not  pure  geometric.  When  pressure  gradient 
and  density  of  matter  decrease  some  coalescence  of  meshes  takes 
place.  . 

Using  moving  grids  demonstra1i,es  a  substantial  economy  of 
reserved  computer  memory  In  comparison  with  a  fixed  eulerlan 
grid.  One  can  calculate  that  about  30  000  -  40  000  meshes  Is 
required  to  perform  computation  of  this  problem  using  eulerlan 
grid  with  the  same  accuracy. 

4.  Hypervelocity  impact  of  iron  projectile  on  a  target 
with  a  conical  hole 

This  example  demonstrates  possibilities  of  moving  grid 
adaptations  to  region's  form  In  the  case  of  strong  deformations. 
In  the  beginning  of  Interaction  shock  pressure  of  about  0.5  Mbar 
Is  generated  In  PMMA  disk,  while  pressures  In  Iron  target  reach 
about  2  Mbar.  Then  two  jets  of  liquid  PMMA  are  formed  moving 
along  the  walls  of  conical  channel  (fig.  16).  Collision  of  the- 
jets  on  the  axis  leads  to  cumulation  effect.  Pressure  behind  the 
projectile  reaches  to  2  Mbar  (fig. 17).  After  Impacting  With  the 
Iron  projectile  these  jets  breach  between  the  projectile  and  the 
walls  of  the  channel  accelerating  projectile  (fig. 18).^ 
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Fig. 9  Grid  and  presspre  levels. 
1-100  Kbar,  2-  200  :;kbar  .  .  .  . 
Time:  710  ns  6'. -i'  '  T  ■  i.i  ;  -  ■ 


Fig. 10  Grid  and  pressure  levels 


1 -too  kbar ,  2-200  kbar. . . 


Z,  CM 


Fig. 11  Grid  and  pressure  levels  (best  accuracy) 
step  100  kbar,  time:  710  ns 


-0.5  0.0  0.5  1.0 

_ _  v  '  •  Z,  CM. 

Fig.  12  Grid  velocity  levelsJ^V10.6/2-^o’37  3-9.9,  4-9.6, 
5-9.2,  6-8.8, 7-8.5,  8-8.1,  9-7.7,  10-7.4  km/s 

Time:  1900  ns 


Fig.  13  Grid  and  velocity  levels  (the  same  as  for  Ilg.i2) 

Time:  1900  ns 
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1- 300  kbar 

2- 600  kbar 

3- 900  kbar 

4- 1200  kbar 

5- 1500  kbar 

6- 1800  kbar 
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Introduction 


In  this  report  we  present  applications  of  previously  described  moving  grid  Godunov  code  to 
several  hypervelocity  impact  problems.  This  method  employs  the  exact  solution  of  Riemann 
problem,  therefore  it  is  well  suited  for  construction  of  "nonhomogeneous"  computational 
algorithms  fitting  flow  singularities.  It  is  clear,  that  to  fit  all  singularities,  which  can  arise  in  the 
flow  is  rather  difficult  algorithmic  problem  esspecially  in  multidimentional  case.  That  is  why  the 
most  of  widespread  computational  methods  of  shock  wave  dynamics  are  based  on 
homogeneous  algorithms  treating  in  the  same  fashion  both  continuos  and  discontinuos  domains 
of  flow. 

To  better  resolve  discontinuities  so  called  high  resolution  shock-capturing  schems  are  used. 
Most  of  them  are  second  order  modifications  of  Godunov  scheme  [1],  which  use  approximate 
Riemann  problem  solvers  [2-5].  One  should  mention  in  this  context  so  called  ENO  schemes 
(essentially  non-oscillatory  schemes)  which  allows  one  to  increase  the  order  of  accuracy  up  to 
fourth  or  even  more  [6-7].  Of  course,  the  higher  is  the  order  of  accuracy,  the  more 
complicated  is  the  computational  algorithm.  For  example,  as  it  was  found  in  [8],  the  cost  of  the 
third-order  ENO  scheme  is  four  times  that  of  the  second-order  one.  It  may  appear  the  second 
order  schemes  to  be  more  profitable  when  using  refined  mesh  than  the  fourth-order  ones. 

We  have  chosen  another  approach  to  construction  of  high  resolution  shock-wave  codes, 
namely,  Godunov  scheme  on  a  moving  grid.  The  main  principles  of  this  approach  were 
described  in  [9-10]  and  in  the  previous  HEDRC  report.  Note  that  in  the  framework  of  moving 
grid  code  shock  waves,  contacts  between  different  materials  and  matter-vacuum  boundaries, 
can  be  easily  fitted  .  By  that  the  grid  boundaries  coinsided  with  these  discountinuities  are 
moved  in  accordance  with  the  solution  of  Riemann  problem. 

In  the  present  report  we  demonstrate  the  advantages  of  moving  grid  approach  for  solving 
four  different  problems. 

The  first  one  is  an  impact  of  thin  plane  projectile  on  a  thin  screen.  A  great  difference  of 
longituidal  and  transverse  dimentions  leads  to  a  great  number  of  zones  when  using  a 
Lagrangian  grid,  which  covers  all  the  projectile  and  the  target.  Our  method  allows  one  to 


compute  only  those  domains  where  the  flow  is  two-dimentional,  that  gives  a  substantional 
economy  of  computer's  resources. 

The  second  example  is  the  numerical  simulation  of  penetration  of  a  cosmic  body  moving  with 
hypervelocity  into  gaseous  atmosphere  with  exponentionally  increasing  density.  The  typical 
feature  of  this  problem  is  the  presence  of  various  space  scales.  For  examples,  the  size  of  the 
body  is  much  larger  than  the  thickness  of  compressed  gas  layer  at  the  face  side  of  the  body  and 
on  the  other  side  it  is  much  less  than  the  distance  which  this  body  passes  before  it  explodes. 

The  third  problem  is  a  debris  cloud  simulation  produced  by  ball-plate  impacts.  In  this 
problem  the  ball  and  the  plate  have  thicknesses  of  the  same  order,  that  makes  it  difficult  to 
employ  any  Lagrangian  code  because  of  the  great  deformations  of  material.  On  the  other  side, 
it  is  known  that  Eulerian  codes  give  a  great  dispersion  of  mass  due  to  averaging  of  flow  values 
over  the  fixed  zone  intervals.  Application  of  the  moving  grig  code  in  which  the  computational 
grid  moves  with  material  near  the  boudaries  and  has  a  minimal  advection  error  in  the  central 
grid  region  allows  one  to  overcome  these  diffuculties. 

The  last  example  is  a  numerical  simulation  of  detonation  of  charges  of  finite  diameters.  It  is 
well  known,  that  a  steady  detonation  can  exist  only  if  the  size  of  the  explosive  is  greater  than 
some  minimum  value.  To  find  this  size  for  a  charge  of  some  specified  geometry  numerous 
numerical  computations  are  required.  When  this  size  is  close  to  the  critical  one  the  time  which 
is  nessesary  for  setting  up  a  quasi  steady  detonation  regime  is  very  long.  It  means  the  we  must 
be  able  to  track  the  propagation  of  the  front  over  a  long  distance.  As  the  region  of 
computation  increases  with  time,  a  grid  refinement  is  required  to  maintain  some  specified  level 
of  computation  accuracy. 

Numerical  simulation  of  impact  of  projectiles  with  a  great  D/L  ratio 

Impactors  with  a  great  D/L  ration  (where  D  is  the  diameter  and  L  is  the  thickness)  are  used 
in  shock-wave  experiments  in  order  to  simpMfy  their  interpretation  reducing  the  last  one  to  ID 
case.  Nevertheless,  the  flow  becomes  two  dimentional  when  the  rarefaction  waves 
fi'om  the  periphery  reach  the  center  of  the  projectile.  To  account  2D  effects  at  the  late  stages  of 
interaction  2D  simulation  is  nessesary.  In  the  present  report  we  demonstrate  one  example  of 


application  of  moving  grid  to  effective  solution  of  such  a  problem.  In  fig.  l(a-d)  numerical  grid 
and  levels  of  constant  density  are  presented  for  high  velocity  impact  of  1mm  thick  aluminum 
disk  on  2mm  thick  aluminum  foil  at  8  km/s  velocity.  The  diameter  of  the  disk  is  4  cm 
In  the  beginning  of  computation  we  specify  two  small  numerical  regions.  One  of  them 
represents  the  edge  of  the  disk  and  the  other  corresponds  to  the  fragment  of  the  foil,  which  is 
in  contact  with  this  edge.  Then  we  move  the  boundary  segment  separating  the  ID  and  2D  flow 
regions  with  characteristics  velocity  towards  to  the  axis  of  simmetry.  The  other  fragments  of 
the  boundary  are:  free  boundaries  (back  sides  of  the  disk  and  the  foil),  contact  between  the 
disk  and  the  foil  and  the  front  of  the  shock  wave  runing  away  from  the  center.  Rarefaction 
waves  coming  from  the  back  surfaces  of  the  projectile  and  the  target  lead  to  an  expantion  of 
shocked  material  and  to  generation  of  tensile  stresses.  That  is  why  the  contact  between  the  disk 
and  the  target  disappears  .  The  type  of  the  boundary  segments  is  changed  automatically  from 
contact  boundary  to  free  one  (fig.  lb)  For  the  later  times  the  size  of  the  flow  under 
consideration  is  much  more  than  initial  one  (see  fig.  Id).  The  straight  line  fragment  of  the 
boundary  separates  the  flow  region  where  the  transversal  rarefaction  waves  are  important. 

In  addition  to  this  example  one  should  note,  that  there  is  a  lot  of  problems  containing 
different  time  and  space  scales  which  could  be  effectively  solved  using  moving  grid  code,  for 
example,  investigation  of  thin  foil  acceleration  with  laser  and  ion  beams,  perforation  of 
multilayered  spaced  shields  and  so  on. 

Computation  of  large  asteroid  penetration  in  Jupiter's  atmosphere 

The  problems  of  asteroid  hazard  are  of  great  interest,  because  the  number  of  ateroids  moving 
in  vicinity  of  Earth  is  increasing  with  time.  An  impact  of  pereodic  comet  Shoemaker-Levy  9 
on  Jupiter  took  place  recently  (in  July  1994).  The  consequences  of  this  impact  were  observed 
both  from  the  Earth  and  space  sattelits.  In  this  report  we  demonstrate  an  example  of  numerical 
simulation  of  the  penetration  of  one  comet  fragment  of  1  km  diameter  into  Jupiter's 
atmosphere.  The  asteroid  was  approximated  by  a  spherical  incompressible  body  moving  with  a 
velocity  of  60  km/s.  Due  to  its  high  velocity  it  passses  a  distance  much  more  than  its  diameter 
before  it  explodes.  As  the  density  of  the  atmosphere  grows  with  distance  the  flow  is  unsteady. 


Extremly  high  velocity  of  impact  leads  to  heating  and  ionization  of  a  gas  flow  after  shock 
compression.  Since  the  flow  behind  the  shock  front  is  responsible  for  ionization  and  luminosity 
of  gas,  which  can  be  observed  an  accurate  resolution  in  this  region  is  required.  It  is  clear,  that 
application  of  Eulerian  codes  to  this  problem  will  lead  to  a  great  diffusion  of  shock  front  and 
contact  interfaces.  Eulerian  computations  also  require  a  great  number  of  computational  zones 
to  cover  the  regions  in  which  the  asteroid  moves.  That  is  why  Eulerian  computations  of  this 
problem  were  performed  in  'reverse  ballistic'  sense  using  an  atmosphere  moving  towards  an 
initially  stationary  fragment  [11]. 

To  avoid  computations  in  undisturbed  domains  we  surrounded  the  asteroid  region  by  a  thin 
region  of  the  gas,  whose  boundary  moves  with  the  characteristics  velocity.  The  size  of  the  gas 
region  grows  in  time  but  we  can  exclude  from  consideration  those  parts  of  the  gas  region, 
which  are  far  from  the  asteroid  and  are  not  of  our  interest.  The  results  of  computation  are 
shown  in  fig.2-4.  In  the  beginning  (fig.2)  of  penetration  the  gas  detaches  from  the  back  side  of 
asteroid  and  a  gas  jet  appears  moving  in  the  opposide  direction.  Then  we  cut  off  the  part  of  gas 
region  corresponding  to  this  jet  (fig.3).  Fig.4  presents  the  results  of  simulation  for  the  time 
when  the  asteroid  has  passed  about  80  km  in  atmosphere.  The  pressures  at  the  frontal  surface 
of  the  asteroid  become  about  0.3-0. 5  kbar,  that  are  close  to  the  strength  of  material.  Further 
computations  require  to  account  for  the  fracture  of  matterial.  In  accordance  with  [12], 
fragmentation  of  the  asteroid  leads  to  a  growth  of  the  heat  flux  inside  it,  because  of  signifficant 
increase  of  the  effective  surface  area.  This  gives  rise  to  the  rapid  transformation  of  the  asteroid 
material  from  the  condensed  to  gaseous  state.  If  we  suppose,  that  the  fracture  and 
fragmentation  of  asteroid  is  due  to  the  stress  gradient  at  its  surface,  we  obtain  the  distance  of 
about  100  km,  where  the  explosion  takes  place.  This  agrees  with  recent  observations. 


Numerical  simulation  of  propagation  of  debris  cloud  produced  by 
ball-plate  impact. 

The  results  of  a  computation  of  a  20  g  spherical  lead  projectile  striking  a  lead  plate  is 


presented  in  Fig.5(a,b).  The  impact  velocity  for  this  problem  is  6.6  km/s.  Fig.6  represents  the 
experimental  X-ray  photograph  of  the  same  problem  at  the  time  moment  30  |a  m  presented  in 
[13],  Simulation  results  are  in  accordance  with  the  experiment.  The  number  of  zones  used 
grows  during  the  computation  to  ensure  some  specified  level  of  accuracy.  The  criteria  for  grid 
refinement  are  not  purely  geometric.  When  the  pressure  gradient  and  the  density  decrease 
some  coalescence  of  zones  takes  place. 

Using  moving  grids  demonstrates  a  substantial  economy  of  required  computer  memory  in 
comparison  with  a  fixed  Eulerian  grid.  One  can  estimate  that  about  10  000  -  20  000  meshes 
is  required  to  perform  computation  of  this  problem  using  Eulerian  grid  with  the  same 
accuracy. 

Simulation  of  detonation  in  high  explosive  charges  of  finite  diameter. 

The  failure  detonation  problem  is  the  problem  of  minimal  explosive  charge  diameter  when 
a  self-sustained  detonation  can  exist.  We  present  below  an  example  of  determination  of  this 
diameter  using  numerical  simulation.  The  flow  is  described  by  Euler  equations.  The  only 
difference  with  hypervelocity  impact  problem  is  the  appearance  in  the  energy  equation  of  a 
source  term  which  is  responsible  for  heat  release.  To  check  whether  the  detonation  becomes 
steady  or  not  it  is  necessary  to  calculate  the  evolution  of  detonation  process  along  the  distance 
of  at  least  several  (perhaps  tens)  diameters.  As  the  downstream  flow  behind  the  shock  is 
determined  by  the  chemical  reaction  kinetics  and,  in  particular,  by  the  value  of  energy 
release,  the  shock  must  be  calculated  accurately. 

The  following  computational  strategy  is  chosen  in  accordance  with  the  problems  listed  above. 
The  computational  region  represents  a  curvelinear  quadrangle.  The  sides  of  the  quadrangle  are: 
the  leading  shock,  the  segment  of  contact  between  detonation  products  and  vacuum,  the 
segment  of  axis  of  symmetry,  and  the  forth  one  is  a  segment  of  straight  line  between  the  free 
surface  and  the  axis  of  symmetry,  which  is  perpendicular  to  the  axis  of  symmetry.  The  velocity 
of  this  side  is  assumed  to  be  directed  along  the  axis  and  to  be  equal  to  max  (u+a,D),  where  D 
is  the  shock  velocity  at  the  axis  of  symmetry  and  u+a  is  the  velocity  of  characteristic 


surface  moving  along  the  axis  and  calculated  using  parameters  of  the  cells  adjacent  to  this 
side.  It  means  that  the  boundary  moves  relative  to  the  matter  ahead  of  it  at  least  with  the  sonic 
velocity.  Therefore  the  flow  parameters,  which  determine  the  fluxes  of  mass,  momentum 
and  energy  throughout  this  side  are  taken  from  internal  cells.  A  pressure  of  200  kbar,  normal 
density  and  zero  velocity  were  taken  as  an  initial  data  for  calculations.  It  was  found  from 
simulations,  that  the  critical  diameter  is  somewhere  between  2.4  and  3  mm.  Results  of 
simulations  corresponding  to  these  two  cases  are  shown  in  Fig.7(a,b)  In  the  case  of  2.4  mm 
diameter  detonation  decays.  Calculations  for  this  diameter  were  performed  with  different 
initial  pressures.  In  all  the  cases  a  steady  detonation  was  not  obtained.  This  fact  allows  to 
draw  a  conclusion  that  the  failure  diameter  of  TNT  lies  within  the  interval  2.4  and  3  mm. 
This  also  agrees  with  the  experiment  [14]. 

CONCLUSIONS 

We  have  demonstrated  the  robotness  and  flexibility  of  developed  2D  code 
for  hypervelocity  problem  computations.  The  main  advantage  of  this  code 
in  comparison  with  Eulerian  high-order  accuracy  hydrocodes  is  an  accurate 
treatment  of  multimaterial  interfaces.  On  the  other  side,  arbitrary 
Lagrangian-Eulerian  methods  (ALE),  which  utilize  lagrangian  motion  of 
interfaces  and  permit  an  arbitrary  mesh  motion  inside  the  computation 
region  are  not  always  adapted  to  dynamically  evolving  interface  shape.  It  is 
the  main  reason  to  develope  ALE  codes  on  unstructured  grids  [15]. 

We  overcome  these  difficulties  with  the  help  of  decomposition  of  numerical 
region  onto  subregions.  This  decomposition  is  provided  automatically  during 
the  computation.  By  that  one  can  govern  this  process  excluding  from 
computations  the  regions  which  are  not  of  our  interest. 
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Fig.l(a-d)  Numerical  grid  at  different  tii 

a) Time  is  51  ns,  density  level  step  is  0.5 

b) Time  is  1000  ns 

c) Time  is  2010  ns 
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Flg-2.  Numerical  grid  and  pressure  levels.  Pressure  level  step  is  0.001  kbar, 
maximum  pressure  is  0.01  kbar.  Time  is  55  (is. 
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Fig.3  Numerical  grid  and  pressure  levels.  Pressure  level  step  is  0.001  kbar, 
maximum  pressure  is  0.014  kbar.  Time  is  80  /is.  Initial  conflguration  is  shown  on 
the  right  side. 
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Fig.4.  Numerical  grid  and  pressure  levels.  Pressure  level  step  is  0.06  kbar, 
maximum  pressure  is  0.6  kbar.  Time  is  1.3  s. 
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Fig.5  Numerical  grid  and  levels  of  constant  density. 

a) Time  Is  15  jxs,  level  step  Is  0.05  g/cc 

b) Time  is  30  /ts,  level  step  is  0.01  g/cc 


Fig.  6  Experimental  X-ray  photograph  at  the  time  30  /ts 
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Introduction 


One  of  the  most  important  problem  of  numerical  solution  of  nonlinear 
hyperbolic  conservation  equation  is  to  obtain  high  accuracy  of  the  solution  in 
both  discontinous  and  continuos  regions.  In  recent  years,  a  lot  of  so  called 
high  resolution  schemes  have  been  proposed,  for  example,  second  order  TVD 
scheme  of  Roe  [1,2]  and  Sweby  [3];  upwind  TVD  of  Marten  [4]  and  Yee  et 
al.  [5];  symmetric  TVD  of  Chakravarthy  and  Osher  [6];  ENO  scheme  of  Marten 
et  al.  [7,8],  FCT  scheme  of  Zalesak  [9]  McDonald  and  Ambrosiano  [10]; 
MUSCL  schemes  of  van  Leer[11],  Goodman  and  LeVeque  [12],  Davis  [13], 
Colella  [14],  and  PPM  of  Colella  and  Woodward  [15].  The  most  of  schemes 
mentioned  above  were  originally  derived  for  the  case  of  perfect  gas  and  for 
ID  geometry  (or  2D  for  fixed  Eulerian  grids).  For  numerical  simulations  of 
hypervelocity  impact  problems,  when  realistic  physical  models  are  used,  we 
face  with  additional  problems.  This  is,  for  example,  a  variety  of  matter 
properties  in  different  phase  states  (solid,  liquid,  gas,  plasma)  ,  when  the 
compressibility  of  matter  differs  up  to  several  orders  of  magnitude.  This 
difference  often  gives  some  unphysical  values  of  variables  appearing  due  to 
errors  of  numerical  approximation  (negative  internal  energy  or  density),  which 
lead  to  a  failing  of  numerical  scheme.  As  it  was  found  in  [16],  some  of 
Godunov-type  second  order  methods  based  on  approximate  Riemann  solver 
are  not  positively  conservative  or  positively  conservative  under  certain 
conditions.  They  can  produce  non-physical  states  with  negative  density  or 
internal  energy  in  rarefaction  waves.  We  based  our  code  on  the  Godunov's 
method  of  numerical  integration  of  Euler  equations  on  a  time  dependent  (so 
called  moving)  grid  [17].  As  it  was  shown  in  [16],  this  scheme  is  positively 
conservative  for  at  least  in  the  case  of  ideal  gas  EOS. 

Let  us  remind  the  main  principles  of  this  code.  In  the  framework  of  this 
approach  only  the  boundaries  of  numerical  grid  are  moved  in  Lagrangian 
fashion.  The  location  of  interior  grid  nodes  are  determined  during  grid 
generation  procedure  by  using  of  coordinates  of  boundary  nodes.  As  the 
material  flows  through  the  grid  cells  there  is  some  diffusion  error  as  well  as  in 
any  Eulerian  method.  Nevertheless  there  exist  several  possibilities  to  reduce 
such  errors  by  generating  "near-Lagrangian"  grids.  In  the  present  version  of 
2D  hydrocode  we  utilize  decomposition  of  numerical  region  onto  subregions 
with  Lagrangian  boundaries,  in  which  grids  are  generated  independently 
using  conformal  mapping  procedure.  This  decomposition  is  provided 
automatically  during  the  computations.  When  every  subregion  consist  of  only 


one  cell.  This  approach  resembles  Lagrangian  method  on  unstructured 
quadrilateral  mesh.  On  the  other  side,  such  a  grid  enables  all  advantages  of  a 
regular  grid.  It  is  known  that  the  numerical  approximation  of  solution  on  an 
orthogonal  grid  is  always  better  than  that  one  on  an  arbitrary  disturbed  grid 
[17].  The  grid  generation  procedure  can  be  governed  by  moving  the 
boundary  nodes  along  the  boundary.  Using  such  a  redistribution  of  the 
boundary  nodes  orthogonal  grids  are  generated  in  subregions.  The  other 
advantage  of  orthogonal  grids  is  a  possibility  to  simplify  a  second  order 
extension  of  Godunov's  scheme. 

In  the  present  report  we  consider  a  second  order  extension  of  moving  grid 
code  and  compare  results  of  calculations  of  hypervelocity  impact  problems 
with  shock-fitting  calculations.  The  strategy  of  shock  fitting  grids  was 
described  in  the  previous  report.  Unfortunately,  the  application  of  this 
approach  is  not  possible  inside  of  numerical  subregions.  It  means,  that  the 
description  of  shock  wave  propagation  and  interaction  inside  of  disturbed 
subregion  can  not  be  done  using  shock  fitting  grids.  That  is  the  reason  why 
the  second  order  extension  of  Godunov's  scheme  is  very  important. 

3.NUMERICAL  ALGORITHM 

In  this  report  we  describe  algorithm  of  computations  only  briefly.  More 
detailed  description  one  can  find,  for  example,  in  [18],  as  well  as  in  the  first 
report.  As  we  have  already  mentioned,  the  numerical  region  is  divided  into 
several  subregions  with  Lagrangian  boundaries  during  the  computations.  Let 
us  consider  the  main  steps  of  computational  algorithm  for  one  subregion. 

The  first  step  of  computations  is  the  displacement  of  the  boundary  of 
subregion.  After  shifting  the  boundary  to  a  new  position  some  segments  of 
the  boundary  can  change  their  type  due  to  interaction  with  boundaries  of  the 
other  subregions.  The  boundary  type  determines  the  boundary  condition 
(what  is  necessary  for  Riemann  problem  initial  data)  and  the  law  of  motion  of 
the  boundary.  For  example,  a  shock-front  boundary  moves  according  to 
Hugein's  principle,  rigid  wall  boundary  does  not  move  and  so  on. 

The  second  step  is  an  iterative  procedure  of  orthogonal  grid  generation 
inside  of  subregion.  This  procedure  is  based  on  the  distribution  of  the 
boundary  nodes.  This  distribution  is  performed  several  times  until  the  grid 
generated  with  conformal  mapping  becomes  orthogonal. 

The  third  step  is  a  solution  of  Riemann  problems  for  inner  zones  and 
calculation  of  fluxes  between  the  neighboring  zones.  One  should  remark  here. 


that  for  an  arbitrary  equation  of  state  the  Riemann  problem  can  be  solved 
only  numerically.  Nevertheless,  we  employ  an  exact  solver  of  this  problem 
only  in  vicinity  of  flow  discontinuity.  The  most  of  Riemann  problem 
computations  are  done  either  in  isentropic  or  in  acoustic  approximation. 

Second  order  extension  of  Godunov's  scheme  can  be  obtained  if  we 
assume  a  piecewise  linear  distribution  of  flow  parameters  inside  of  the  grid 
cells.  To  conserve  the  monotonicity  property  of  Godunov's  scheme  we  use 
the  "minimum  derivative  principle"  proposed  in  [19].  The  main  idea  of  this 
principle  is  to  choose  the  minimum  possible  derivative  when  interpolating  the 
values  from  the  zone  center  to  the  boundary  with  the  neighboring  zone, 
where  the  Riemann  problem  must  be  solved.  If  the  grid  is  orthogonal, 
derivatives  only  in  one  direction  (either  along  the  grid  rows,  or  along  the  grid 
columns)  are  important  when  interpolating  cell-centered  values  to  the 
boundaries  of  the  cells.  This  simplifies  the  realization  of  "minimum  derivative 
principle" . 

4. EXAMPLES  OF  CALCULATIONS 

In  the  present  paper  we  consider  an  application  of  our  2D  code  to 
hypervelocity  impact  problem.  One  of  the  most  difficult  task  for  any  2D 
Lagrangian  code  is  to  compute  hypervelocity  penetration  of  projectile  into  a 
thick  target.  We  consider  two  different  approaches  to  numerical  simulation  of 
hypervelocity  impact  of  iron  projectile  moving  with  10  km/s  velocity  on  a 
thick  aluminum  target  . 

First  one  concludes  in  a  usage  of  "near-Eulerian"  grid,  which  covers  the 
region  of  the  target  where  the  shock  wave  propagates. 

The  second  one  is  fitting  of  this  shock  wave  by  moving  the  boundary  of  the 
grid  with  shock  front  velocities  obtained  by  Riemann  solver.  The  second 
method  in  fact,  is  more  accurate  to  resolve  the  shock  front  and  is  more 
economic,  because  it  does  not  require  a  grid  in  undisturbed  regions.  On  the 
other  hand  the  shape  of  numerical  region  is  more  complicated  and,  as  a 
consequence  of  that,  the  grid  generated  is  not  so  orthogonal  as  in  the  first 
case.  The  results  of  computations  corresponding  to  these  cases  together  with 
numerical  grids  are  shown  in  Fig.1.  There  were  two  subregions  used  in  the 
case  of  shock-fitting  grid.  One  of  them  corresponded  to  the  projectile,  the 
other  one  covered  the  disturbed  region  in  the  target.  In  the  other  case  (near- 
Eulerian  grid)  three  subregions  were  used  in  computations. 


Fig.1  Numerical  grid  and  pressure  levels  for  hypervelocity  impact  calculations 
using  shock  fitting  grid  (above  the  axis)  and  "near-Eulerian"  grid  (below  the 
axis). 

It  is  known,  that  first  order  Godunov's  scheme  is  a  very  dissipative  one. 
Our  experience  shows,  that  this  scheme  smears  out  weak  discontinuities  in  the 
case  of  realistic  EOS  much  stronger  than  in  the  case  of  perfect  gas. 
Moreover,  this  dissipation  grows  in  time.  Results  of  computation  of 
hypervelocity  impact  of  iron  projectile  on  aluminum  target  shown  in  fig. 2-3 
confirms  this  fact.  In  contrast  to  the  first  order  scheme,  second  order 
calculations  exibit  nearly  constant  width  of  shock  front  smoothing. 


Fig. 2  Grid  boundaries  and  levels  of  constant  pressure  for  the  first  (below  the 
axis)  and  second  order  accuracy  calculations  of  hypervelocity  impact  of  iron 
projectile  (on  the  right)  on  an  aluminum  target. 


Fig. 3  Grid  boundaries  and  levels  of  constant  pressure  for  the  first  (below  the 
axis)  and  second  order  accuracy  calculations  of  hypervelocity  impact  of  iron 
projectile  (on  the  right)  on  an  aluminum  target. 


As  it  has  been  shown  in  [20]  for  the  case  of  steady  gas  flow,  the  second 
order  extension  of  Godunov's  scheme  is  not  so  sensitive  to  distortions  of 
numerical  grid  as  it's  first  order  counterpart.  We  demonstrate  below  (Fig. 4-5) 
the  results  of  simulation  of  hypervelocity  impact  obtained  using  the  first  and 
the  second  order  Godunov  scheme  both  on  the  shock  fitting  grid  and  on  the 
near-Eulerian  grid.  Our  results  generally  confirm  the  conclusions  obtained  in 
[20].  For  the  second  order  scheme  (Fig. 5)  the  difference  between  the  results 
obtained  on  the  different  grids  is  negligible,  while  for  the  first  order  scheme 
there  is  a  big  difference  in  the  amplitude  of  the  shock  wave  calculated  on  the 
shock  fitting  and  on  near-Eulerian  grid.  It  is  naturally,  that  the  shock  front 
dissipation  is  much  pronounced  in  the  case  of  the  first  order  scheme  than  in 
the  case  of  the  second  order  one. 


Fig. 4  Grid  boundaries  and  pressure  levels  (in  kbars)  for  shock  fitting  grid  and 
"near-Eulerian"  grid.  First  order  scheme. 


Second  order  scheme  tjjne=7749  ns 


Fig. 5  Grid  boundaries  and  pressure  levels  (in  kbars)  for  shock  fitting  grid  and 
near-Eulerian  grid.  Second  order  scheme. 


CONCLUSIONS 

We  have  considered  an  application  of  Godunov's  scheme  on  the  moving 
grids  and  it's  second  order  extension  (Kolgan  scheme)  on  the  example  of 
hypervelocity  impact  calculation.  We  have  shown  that  this  scheme  coupled 
with  the  algorithm  of  decomposition  of  numerical  region  onto  subregions 
shows  a  good  robustness  and  flexibility  for  multimaterial  problems  of 
computational  fluid  dynamics  with  large  distortions.  The  other  advantage  of 
decomposition  of  numerical  region  of  complicated  shape  is  a  possibility  to 
generate  orthogonal  grids  in  subregions 

We  have  developed  second  order  Godunov's  scheme  applying  Kolgan's 
algorithm  of  "minimal  derivative".  It  has  been  found  that  Kolgan's  scheme  is 
not  so  sensitive  to  the  form  of  the  moving  grid  as  original  Godunov's  scheme. 
The  results  of  second  order  calculations  are  much  closer  to  more  accurate 
shock-fitting  calculations. 
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Brief  description  of  the  code 


The  language  of  the  code  is  standard  F77  FORTRAN.  So  it  can  be  easily 
installed  on  any  computer  with  f77  compiler.  The  code  consists  of  the 
following  program  modules: 


Module:  Function: 


STARTT  -  input  of  initial  data  and  transformation 

of  these  data  into  file  NAM.tmp  for 
continuation  of  calcidations  by  module  BIG2T 

BIG2T-  carries  out  calculations,  renews  data  for 

continuation  of  calculations  (file  NAM.tmp) 
and  records  results  of  the  calculation 


REFT-  modification  of  the  task  (possibihty  to 

add  or  delete  subregions  ).  The  changed  data 
are  rewritten  into  the  file  NAM.tmp 

UNFTOF-  transformation  of  NAM.tmp  from  binary  file  to 
ASCII  file  and  back  (it  may  be  useful  for 
transmission  of  calculation  results  from  one 
computer  to  other  one  with  different  binary 
data  presentation) 


This  version  of  2D  hydrocode  operates  with  tabulated  equations  of  state 
(TABEOS  format)  as  well  as  with  perfect  gas  analytical  EOS.  The  code 
loads  four  EOS  tables  :ENER*.TAB,PRES*.TAB,SND*.TAB, 
PHYS*.TAB,  DATA*.MAN  (*  -  from  1  to  4),  which  are  energy,  pressure, 
sound  velocity,  index  of  state  and  information  about  table  hmits  for  1- 
aluminum,  2-Nickel,  3-  tantalum,  4-  tungsten  correspondingly.  If  it  is 
required  to  include  another  TABEOS  tables  in  computations  one  need  to 
change  the  normal  densities  for  substances  1-4  in  subroutine  readtab 
(urst.for)  in  accordance  with  the  tables  used. 


Installation  of  the  code 

To  create  object  modules 

1)  copy  aU  files  *.FOR  and  *.INC  to 
one  directory  on  a  hard  disk 

REMARK:  all  *. INC  files  must  be  with  capital  characters 


2)  compile  aU  *.FOR  files 


*  before  compilation  probably  may  be  necessary: 

A)  to  set  function  which  determines  CPU  time  in  seconds  in 
SUBROUTINE  ETIME  in  the  file  PC.FOR.  (to  decomment 
correspondent  string  may  be  sufficient) 

B)  to  change  parameters  in  the  file  BIGPAR.INC 
(dimensions  of  arrays) 

3)  create  loading  module  STARTT,  using  object  files: 

startt.o,  blocka.o,  adi.o,  arct.o,  big.o,  contw.o,  univer.o,  datO.o, 
dattst.o,  disto.o,  jup.o,  pc.o,  shr.o,  spect.o,  urst.o 

4)  create  loading  module  BIG2T,  using  object  files: 

basica.o,  blocka.o,  mainxt.o,  spect.o,  tnew2a.o,  urst.o,  adi.o,  arct.o 
big.o,  buil.o,  builO.o,  buill.o,  build.o,  contw.o,  datO.o,  del.o,  disto.o 
god2.o,  gmew.o,  jup.o,  pc.o,  ref.o,  shr.o,  split.o,  split l.o,  univer.o 

5)  create  loading  module  REPT,  using  object  files: 

rep.o,  blocka.o,  spect.o,  urst.o,  adi.o,  arct.o,  big.o,  contw.o,  datO.o 
del.o,  jup.o,  pc.o,  ref.o,  shr.o,  split.o,  sphtl.o,  univer.o 

*Note  that  modules  mentioned  above  (STARTT,  BIG2T,  REPT)  can  be 
created  using  makefile.  To  generate  first  order  accuracy  code  (BIGIT), 
one  should  replace  tnew2a.*  by  tnewla.*  in  the  makefile. 

Generation  of  a  new  initial  data  file 

1)  Chose  one  of  the  existing  initial  data  files  from  the 
file  fist  below 
TST2.DAT 
VRB2.DAT 
VRB3.DAT 


2)  Copy  the  chosen  file  into  file  NAM.dat,  where  NAM  is 
any  name  consisting  of  1-4  symbols 

(the  names  of  aU  the  files  created  during  the  computation 
of  this  problem  will  include  NAM  as  a  root  word) 

3)  Correct  contents  of  the  file  NAM.dat  with  text  editor 
using  description  of  initial  data  file  presented 

below 


Description  of  the  structure  of  initial  data  file  NAM.dat 


The  initial  data  file  consists  of  two  parts: 

First  part  contains  control  information.This  part  has  the  same  format  for  all 
files  of  this  type.  Let  us  consider  an  example  shown  below,  where  the 
fist  part  of  initial  data  file  TST.DAT  is  presented. 

AL/AL  U=10.1KM/C 

JSPR-LOS-LST-IRD-ISP-IRA-ICM-AIR 
0  1  1  2  1  00  1  0 
NBL-MI-ID1-ID2-IK-IGR-ILENPARIVSP 
9900  1  10  10  300  10  -10  6  0 

- CMM— : - H3 - ; - HDEL— : - T - STABL— : 

-l.OOOOOE-03  0.00075E-00  O.OlOOOE+00  l.lOOOOE+05  0.50000E+00 


The  first  string  is  comment.  This  is  usually  short 
information  about  the  problem.  (Here  this  is  10.1  km/s  impact 
of  aluminum  on  aluminum).  The  next  strings  consist  of 
title_strings  and  contents_strings.  The  title_strings  show  the 
names  of  variables.  The  contents_strings  are  values  of  these 
variables.  The  assignment  of  the  variables  is  shown  in  the  table 

TABLE 


variables:  meaning: 

ispr  controls  rebound  of  different  boundary  segments 

0-  calculations  without  rebound 

1-  rebound  for  "contact"  type  boundaries  only 

2-  rebound  for  "contacts  and  rigid  walls"  type 
boundaries 

3-  rebound  for  all  boundary  types 

los  1-  cylhndrical  symmetry 

0-  plane  symmetry 

1st  1-  there  are  rigid  walls 

0-  rigid  walls  are  absent 

ird  0-  to  ignore  isp 

1-  to  continue  calculations  from  the  restart  data 

2-  to  correct  in  dialog  regime  the  number  of  time  steps 
after  which  the  sphtting  of  numerical  regions  into 
subregions  wiU  be  switched  on 

isp  1-  to  change  in  dialog  regime  the  type  of "  contact " 
between  subregions 

0-  to  do  nothing 


ira  energy  source  (for  impact  problems  ira  is  equal  to  0) 

icm  1-  automatical  choice  of  time  step 

0 - time  step  =  abs(cmm)  *  stabl 

(cmm,stabl  see  below.) 

nbl  is  the  number  of  time  steps  to  do  before  copying  of  all 
flow  parameters  to  output  file. 

mi  1-  to  search  the  optimum  location  of  4  angle  tops  of 
rectangular  subregion  along  the  boundary 

idl,id2  are  the  numbers  of  time  steps  to  do  before  the 
data  record  for  restart(idl-  before  the  first  record) 

ok  is  the  maximum  number  of  time  steps  to  do. 

igr  is  the  number  of  time  steps  to  do  before  printing  packed 
information  about  aU  subregion's  boundaries. 

ilen  is  the  number  of  time  steps  to  do  before  the  next  graphic 
data  record. 

If  ilen  less  0  then  data  for  the  graphics  will  be 
written  in  time  step=abs(ilen*cmm)  in  nanoseconds. 

(cmm  parameter  described  below) 

par  is  the  number  of  flow  parameters  which  must  be  written  for 
graphics  data,  (from  5  to  7  parameters) 

ivsp  if  ivsp  is  not  equal  to  0  then  ivsp  is  input  parameter  for 
subroutine  AUTOSP,  which  changes  types  of "  contact " 
between  subregions  automatically  . 

This  automatical  regimes  must  be  programmed  by  user  for  the 
concrete  task  in  subroutine  AUTOSP. 
cmm  is  parameter  which  determines  time  step  (stt) 
stt  =  cmm  *  stabl  if  cmm>0. 

If  cmm<0.  then  cmm  will  be  read  from  restart  or 
cmm=  abs(cmm)  in  case  of  initial  data. 

h3  space  precision.(H3  equals  to  0.1  of  cell  size) 

H3  is  widely  used  in  the  code  as  character  space  size. 

hdel 

The  boundary  section  with  length  <  hdel*  (cell  size)  will  be 
combined  with  neighboring  one. 

t  time  cpu  Umit  in  seconds 


stabl  stability  coefficient  stt=  cmm  *  stabl 


The  second  part  of  the  initial  data  file 
includes  an  information  about  geometric  configurations  of  the 
problem  to  be  solved.  The  information  may  be  prepared 
using  so  called  UNIVERSAL  INPUT.  In  this  case  the  first  string 
in  the  second  part  of  the  initial  data  file  must  be 
UNIVER . 

One  Universal  input  starting  with  UNIVER.. and  finishing  with 
UNIEND  specifies  only  one  subregion  (numerical  block). 

To  specify  several  subregions  one  must  use  this  block  of 
input  strings  several  times. 

The  boundaries  of  numerical  subregion  are  specified  in 
accordance  with  the  following  rules: 

1)  If  the  problem  contains  the  rigid  walls  and  shock  fronts 
they  must  be  specified  right  after  the  string  UNIVER.... 

For  example,  in  the  case  of  rigid  wall  the  input  looks  like 
this 

UNIVER . 

RIGIDRD 

IWALL,  INW,  ITP,  AW,  BW,  CWjSTRAIGHT  L1NE;F=1NW*(AW*X+BW*Y-CW)=0. 

1  ,  -1  ,  1,  0.  ,  1.,  0. 

-1  ,  -1  ,  1,  0.  ,  1.,  0. 

Here  the  second  string  (after  UNIVER)  starting  as  RIGIDRD  indicates 
to  the  beginning  of  wall  specification,  which  is  finished  when  the  first 
parameter  (IWALL)  equals  to  some  negative  value  (-1  in  this  example).  In 
the  example  presented  above  the  rigid  wall  specified  is  the  straight  Une 
(ITP=1)  and  input  parameters  are  the  coefficients  which  determine  it's 
orientation. 

To  specify  the  front  one  must  put  the  string  starting  as  FRNTREAD 
either  after  the  UNIVER...  or  after  the  rigid  wall  specification.  The  front 
specification  looks  Uke  the  following 

FRNTREAD 

cl  —  read  of  front  information 
ifr,ityp,  |  the_number_of_front,  the_type_of_front 
icont,  1  the_number_of_contour 

a,b,c,  I  the  contour  is  straight  line  a*x+b*y-c=0  !  icont=nomco+l 
ibeg  end,]  l-begin,2-end 
-1 


Here  the  expressions  on  the  right  from  |  are  comments.  The  front 
input  is  finished  when  ifr  equals  to  some  negative  value  (-1  in  this 
example).  The  beginning  and  the  end  of  the  front  move  along  the  stright 


lines  (called  contour).  There  are  not  more  than  two  contours  corresponding 
to  the  same  front.  The  sequence  of  the  numbers  of  the  contours  specified 
must  be  increasing.  Every  contour  corresponds  either  to  the  beginmng  of 
the  front  (ibeg_end=l)  or  to  the  end.  Using  such  an  input  several  fronts 
can  be  specified.  There  are  several  front  types  considered 
(the_t)pe_of_front).  The  main  of  them  are 

1-  shock  front  in  imdisturbed  solid  media 

2- shock  front  in  undisturbed  porous  media 

Specification  of  all  fronts  and  walls  must  be  closed  by  the  string 

END . 

2)  The  next  two  strings 

5.,  0.,  0.,::-USRT-,— VSTR— ,VPOL 
12,::-NUMBER  OF  URS 

give  the  velocity  components  (USRT,  VSTR),  porosity  of  the  subregion 
(VPOL)  and  number  of  material  (EOS).  Then  the  boundary  segments 
must  be  specified.  Boundary  orientation  is  counter-clockwise.  The  boundary 
segment  can  be  straight  fine,  circle  or  some  other  arbitrary  curve  specified 
in  the  table  form.  The  next  string  gives  the  coordinates  of  the  first  point 
of  the  boundary. 

0.,0.,::X0,Y0  -  the  location  of  the  first  point  for  a  boundary  (cm) 

If  the  first  boundary  segment  is  a  straight  line,  then  the  next 
two  strings  are 

STRAIGHT  1 

0.,1.,  XE,YE,::  IS  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

If  the  next  segment  is  a  fragment  of  a  circle  then  must  be 
CIRCLE 

XC,YC,ALF,::  is  the  center  of  the  CIRCLE(CM)^GLE  (0-360  grad) 

The  boundary  segment  can  be  any  arbitrary  curve.  In  this  case  it 
is  specified  using  the  table  input 

TABLE 

XE,YE,::  IS  THE  NODES  COORDINATES 

.......  ,  TABLE  INPUT  UNTIL  INPUT  ERROR  OCCURES 

END  OF  THE  TABLE 

The  specification  of  the  boundary  of  a  subregion  is  closed  by 
the  string 


REGEND 


The  four  angles  of  subregion  must  be  determined  for  the  purpose  of  grid 
generation,  when  the  boimdary  of  the  subregion  is  specified.  The  angles 
input  looks  hke  this 

1___2— 3— 4  ANGLES  LOCATION 

11,12,13,14,;;  to  find  4  angles  of  the  rectangular  region  using 

CORRESPONDENT  BOUNDARY  SEGMENTS  NUMBERS  11,12,13,14 

In  the  example  of  input  shown  below 

the  ANGLES  with  numbers  1,2,4  are  determined  automatically  and  the 
angle  3  is  found  in  the  beginning  of  the  boundary  segment  number  5, 
Initially  all  the  boundary  segments  have  the  type  of  "FREE  BOUNDARY". 
If  one  needs  to  assign  some  other  boundary  type  to  the  boundary  segment 
the  following  strings  are  necessary; 

BOUNDARY  TYPE  INPUT 
Then  depending  on  the  boundary  type  can  be; 

a) in  the  case  of  rigid  wall 

RIGID  ,;CHANGES  "FREE"  TYPE  BY  "RIGID" 

IREG,IRG  ,;  ireg-number  of  boundary  segment 
IRG-NUMBER  OF  RIGID  WALL 

b) in  the  case  of  shock  front 

FRONT  ,;CHANGES  "FREE"  TYPE  BY  "FRONT’ 

IREG,IFR  ,;  ireg-number  of  boundary  segment 
ifr-number  of  front 

c)  in  the  case  of  contact  discontinuity 
CONTACT  ,;cHANGES  -free-  type  by  -contact’ 

IREG,N1,IREG1,;  IREG-NUMBER  OF  BOUNDARY  SEGMENT  FOR  CURRENT  SUBREGION,  Nl- 
NUMBER  OF  SUBREGION  WHICH  IS  IN  CONTACT,  IREGl-NUMBER  OF  BOUNDARY 
SEGMENT  FOR  N1  SUBREGION 

The  input  of  data  for  the  current  subregion  is  closed  by 
UNIEND 
Remark; 

We  have  described  UNIVERSAL  INPUT  in  the  text  above.  If  there  are 
special  program  for  more  simple  input  for  concrete  configuration  then  the 
name  of  this  program  is  the  first  string  in  the  second  part  of  the  input  data. 
For  example; 

DATSTl  ;  -name  of  the  program  which  treats  input  data 


How  to  run  the  code 


1)  Run  the  module  STARTT  in  the  directory  containing  file 
NAM.dat 

2)  Type  NAM  when  the  question  "ENTER  DATA  NAME"  appears 

module  STARTT  will  create  the  files: 

NAM.tmp-  for  following  calculations 
NAMOOO.dat-  first  file  for  graphics 

NAM.cfg-  file  containing  some  parameters  (if  this  file  does  not 
exist) 

An  example  of  CFG  file  is  presented  below 


- -VRB2.CFG - 

0.7654392  ,  PLIM  [bar]  fine  1 

2.7104001E-08,  RLIM  [G/(CM**3)]  line  2 
13.93920  ,  EPSGAS  [M/C]  hne3 

278.7840  ,  EPSCON  [M/C]  hne4 

9.9999998E-03,  ALMIN  hne5 

0.3000000  ,  SINCOND  hne6 

2.0000000E-02,  ALMN  hne7 

9.9999998E-03,  SHISAMI  Hne8 

5.0000001E-02,  SHISAMA  hne9 

0.07000000  ,TET  linelO 

0.000000  ,ADP  hnell 


In  this  file  the  first  two  strings  are  minimal  pressure  and  density,  the  next 
two  strings  are  accuracy  of  calculation  of  wave  velocities  for  the  cases  of  gas 
EOS  and  table  EOS  correspondingly.  Parameters  ALMIN,  SINCOND 
(SINCOND  is  changing  from  0  to  1),  ALMN  are  responsible  for  the 
movement  of  boundary  nodes  on  the  boundary  of  "special"  type.  SHISAMI 
and  SHISAMA  are  minimum  and  maximum  shift  of  the  boundary  nodes 
during  the  procedure  of  boundary  node  distribution  (in  relative  units).  TET 
is  a  parameter  determining  the  convergence  of  grid  generation  (TET>0) 
and  ADP  is  a  level  of  grid  adaptation  to  solution.  In  the  present  version  of 
the  code  adaptation  in  not  included  (ADP=0) 

When  the  STARTT  finishes  its  work  the  following  message  appears 

THE  NUMBER  OF  PICTURE  =0 

3)  Run  the  module  BIG2T  (or  BIGIT)  to  continue  calculations 

4)  Reply  the  questions  appearing  on  the  screen 

To  avoid  the  grid  tangle  in  the  cases  of  large  deformations  to  spht  the 
numerical  regions  into  subregions  is  recommended.  The  most  frequently 
appearing  question  concerning  that  is 


IF  YOU  DON’T  WANT  TO  SPLIT  SUBREGION  * 


ENTER  number  of  steps  without  splitting  analysis 

Enter  0  on  this  question  if  you  want  to  spht  subregion,  which  the  program 
suggests  to  spht.  Then  press  enter  to  continue  computations. 

The  question  ENTER  N,M  ACCORDING  TO  THE  STENCIL  means 
that  you  can  change  the  type  of  contact  between  subregions  N  and  M. 
Enter  -1,1  to  ignor  this  (and  -1,1  for  the  next  question) 

The  program  periodically  creates  the  files  NAMxxx.dat  for  graphics  and 
renews  file  NAM.tmp  for  restart. 

When  the  program  has  done  the  calculations  it  starts  the  dialogue 
from  the  beginning  (asking  "ENTER  DATA  NAME").  To  continue 
the  computations  of  the  same  problem  enter  the  same  name  NAM. 

5)  Analyze  the  results  of  computation  using  graphics 
modules.  For  this  purpose  created  graph  files 
NAM***.dat  (and  kNAM***.dat)  must  be  copied  to  IBM  PC 
Then  they  can  be  depicted  on  the  screen  using  and 
written  in  TIF  file  using  GR.EXE 

(See  the  manual  for  GR.EXE) 

6)  Continue  calculations  (point  3)  if  it  is  necessary 

7)  Change  configuration  of  the  task,if  it  is  necessary 
if  you  want  to  change  configuration  you  must  create 
the  file  dNAM.dat  (see  description  of  the  file) 

and  run  the  program  REFT 

Brief  manual  for  GR.EXE 

Executing  module  GR.EXE  shows  the  results  of  computation  on  the 
screen  of  IBM  PC  (DOS  is  required)  and  writes  the  pictures  shown  on  the 
screen  in  TIF  format.  Besides  the  resulting  files  ROOT_*.DAT  to  be 
shown,  on  more  Ele  K_ROOT.DAT is  necessary  to  determine  the  limits  of 
the  picture.  This  Gle  can  be  created  anytime  by  coping  ateady  existing  Eles 
of  such  type.  After  running  GR.EXE  one  need  to  answer  the  questions  in 
the  dialog. 

1)  Type  the  ROOT  name  of  *. DAT  Files  to  be  plotted 

2) Enter  the  number  of  picture  N  (Ele  ROOT_N.DAT  will  be  plotted) 

3) Choose  the  Eow parameter  to  be  plotted  (1-density,  2-pressure,3- 
porosity,  4, 5-velocities,  6-energy) 

4) Choose  the  type  of  picture  (I-levels  of  constant  values  and  numencal 
grid,  12-levels  of  constant  values  for  two  different  Eow  parameters,  DY- 
cross  sections  along  the  axis  of  symmetry,  DX-  perpendicular  to  the  axis) 


5)  When  the  picture  is  shown  on  the  screen  one  can  plot  over  this  picture 
another  pictures  (if  the  axes  are  not  plotted  yet).  For  this  purpose  enter 
negative  number  -N  (where  N  is  the  number  of  picture  to  be  plotted) 

6)  To  finish  plotting  enter  any  character  (for  example  ')  and  axes  will  be 
shown.  Then  type  enter ,  the  name  of  the  TIF  Hie  where  this  picture  will 
be  saved  and  CTRL-C  to  leave  the  program. 

7)  REMARK:  Wien  using  DY  or  DX  options  (plotting  cross  sections)  use 
CTRL-Z  to  return  back  for  plotting  pictures  after  the  question 

1=  INPUT  Y/X  COORDINATE  FOR  SECTION. 

Example  of  initial  data  files 

Let  us  consider  some  demonstrative  examples  of  initial  data  generation 
for  hypervelocity  impact  problems 

EXAMPLE  1;  initial  data  file  VRB3.DAT 

Impact  of  W  projectile  on  A1  target 
JSPR-LOS-LST-IRD-ISP-IRA-ICM-AIR 
1  1  1  2  1  00  1  0 

NBL-MI-ID1-ID2-IK-IGR-ILENPARIVSP 
9900  1  30  30  300  30-500  6  0 

- CMM— ; - H3 - -HDEL— T - STABL— : 

-l.OOOOOE-03  0.02000E-00  0.06000E+00  l.lOOOOE+05  0.60000E+00 

UNIVER  -NAME  OF  PROGRAM  WHICH  TREATS  INPUT  DATA 
RIGIDRD  ,  ;  READING  AXIS  AND  RIGID  WALLS 

C  IWALL,  INW,  ITP,  AW ,  BW,  CW  ;STRAIGHT  LINE:  F=AW*X+BW*Y-CW=0, 

1  ,  -1  ,  1,  0.  ,  1.,  0. 

-1  ,  -1  ,  1,  0.  ,  1.,  0.  exit, A, A, A, A 

END 

-5.,  0.,  0.,::-USRT-,— VSTR— ,VPOL 

44,::-NUMBER  OF  EOS 

0 . ,  0 . , : :  xo,YO  -  the  location  of  the  first  point  for  the  boundary  (cm) 

STRAIGHT  1 

1.. 0.,  XE,YE,::  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  2 

1. . 4.,  XE.YE,::  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  3 

0.,4.,  XE,YE,::  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  4 

0.,0.,  XE,YE,::  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

REGEND 
c-1-2-3-4 
1,2,3, 4, 

BOUNDARY  TYPE  INPUT 

RIGID  ,;CHANGING  “FREE”  TYPE  BY  "RIGID" 

1,1,;  IREG-NUMBER  OF  BOUNDARY  SEGMENT, IRG-NUMBER  OF  RIGID  WALL 

UNIEND  1 

UNIVER  ,  -NAME  OF  PROGRAM  WHICH  TREATS  INPUT  DATA 


END 

0.,  0.,  0.,::-USRT-,— VSTR— ,VPOL 

41,  ::--NUMBER  OF  EOS 

-8.5,0.,::  xo,yo  -  the  location  of  the  first  point  for  the  boundary  (cm) 

STRAIGHT  1 

0.,0.,  XE,YE,::  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

TABLE  1 

0.,0., 

0.,8.5, 

0.,8.5, 

-8.5,8.5, 

-8. 5, 8. 5, 

-8.5,0., 

-8.5,0., 

END  OF  TABLE 
REGEND 

3— 4 

1,2,„ 

BOUNDARY  TYPE  INPUT 
RIGID 
1,1, 

UNIEND  2 
ENDINPUT 


subregion  Z 


-5.G 


Z.  CM 


-2'.5 


subregion  1 


cf.o 


contact 


Fig.l  Initial  configuration  for  initial  data  file  VRB3.DAT 


EXAMPLE  2:  initial  data  file  VRB2.DAT 
IMPACT  with  shock  fitting 

JSPR-LOS-LST-IRD-ISP-IRA-ICM-AIR 
1  1  1  2  1  00  1  0 


NBL-MI-ID1-ID2-IK-IGR-ILENPARIVSP 
9900  1  30  30  300  30-500  6  0 

- CMM— : - H3— - HDEL— : - T STABL— : 

-l.OOOOOE-03  0.02000E-00  0.06000E+00  l.lOOOOE+05  0.60000E+00 

UNIVER  j  ;  -NAME  OF  PROGRAM  WHICH  TREATS  INPUT  DATA 

FRNTREAD 

Cl  —  read  of  the  front  information 

1,1,  ;c2  ifr,ityp  ,:  the_number_of_front,  the_t5^e_of_front 

1,  ;c3  icont,:  the_number_of_contour 

1.. 0..0.,  ;c4  a,b,c,:  the  contour  is  straight  hne  a*x+b*y-c=0  ! 

icont=nomco+l 

1,  ;c5  ibeg_end,:  1 -begin, 2-end 

-1,  ;c6  end 

RIGIDRD  ,  I  READ  OF  AXIS  AND  RIGID  WALLS 

C  IWALL,  INW,  IIP,  AW  ,  BW,  CW  jSTRAIGHT  LINE:  F=AW*X+BW*Y-CW=0. 

1  ,  -1  ,  1,  0.  ,  1.,  0. 

-1  ,  -1  ,  1,  0.  ,  1.,  0.  exit,A,A,A,A 

END 

-5.,  0.,  0.,::-USRT-,— VSTR— ,VPOL 

44,  ::-NUMBER  OF  EOS 

0.,  0.,  ::  X0,y0  -  THE  LOCATION  OF  THE  FIRST  POINT  FOR  THE  BOUNDARY  (CM) 

STRAIGHT  1 

1  .,0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  2 

1.. 4.,  XE,YE,:;  ARE  COORDINATES  OF  THE  END  OF  THE 
STRAIGHT  SEGMENT  (CM) 

STRAIGHT  3 

0.,4.,  XE,YE,:;  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  4 

0.,0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

REGEND 

1.. .2—3—4  ANGLE  LOCATION 
1, 2,3,4, 

BOUNDARY  TYPE  INPUT 

RIGID  ,;YOU  WILL  CHANGE  "FREE"  TYPE  ON  "RIGID" 

1.1, ;  IREG-NUMBER  OF  BOUNDARY  SEGMENT,IRG-NUMBER  OF  RIGID  WALL 

UNIEND  1 

UNIVER  ,  ;  -NAME  OF  PROGRAM  WHICH  TREATS  INPUT  DATA 

END 

0.,  0.,  0.,::-USRT-,— VSTR— ,VPOL 

41,  ::-NUMBER  OF  EOS 

-0.25,  0.,  ::  X0,Y0  -  the  location  of  the  first  POINT  FOR  THE  BOUNDARY  (CM) 
STRAIGHT  I 

0.,0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  2 

0.,4.1,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

TABLE  3 

0.,4.1,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 


-0.,4.1,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

-0.25,3.0,  XE,YE,;;  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

-0.25,3.0,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) , 

-0.25,0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

-0.25,0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

END  OF  TABLE 
REGEND 

1.. .2—3—4  ANGLE  LOCATION 

1.2.3., 

BOUNDARY  TYPE  INPUT 

RIGID  ,;YOU  WILL  CHANGE  "FREE"  TYPE  ON  "RIGID" 

1.1, ;  IREG-NUMBER  OF  BOUNDARY  SEGMENT, IRG-NUMBER  OF  RIGID  WALL 

FRONT 

3.1, ;  IREG-NUMBER  OF  BOUNDARY  SEGMENT, IRG-NUMBER  OF  FRONT 

UNIEND  2 
ENDINPUT 


Fig.2  Initial  configuration  for  initial  data  file  VRB2.DAT 
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Introduction 


One  of  the  most  important  problem  of  numerical  solution  of  nonlinear 
hyperbolic  conservation  equation  is  to  obtain  high  accuracy  of  the  solution  in 
both  discontinous  and  continuos  regions.  In  recent  years,  a  lot  of  so  called 
high  resolution  schemes  have  been  proposed,  for  example,  second  order  TVD 
scheme  of  Roe  [1,2]  and  Sweby  [3];  upwind  TVD  of  Marten  [4]  and  Yee  et 
al.  [5];  symmetric  TVD  of  Chakravarthy  and  Osher  [6];  ENO  scheme  of  Marten 
et  al.  [7,8],  FCT  scheme  of  Zalesak  [9]  McDonald  and  Ambrosiano  [10]; 
MUSCL  schemes  of  van  Leer[11],  Goodman  and  LeVeque  [12],  Davis  [13], 
Colella  [14],  and  PPM  of  Colella  and  Woodward  [15].  The  most  of  schemes 
mentioned  above  were  originally  derived  for  the  case  of  perfect  gas  and  for 
ID  geometry  (or  2D  for  fixed  Eulerian  grids).  For  numerical  simulations  of 
hypervelocity  impact  problems,  when  realistic  physical  models  are  used,  we 
face  with  additional  problems.  This  is,  for  example,  a  variety  of  matter 
properties  in  different  phase  states  (solid,  liquid,  gas,  plasma)  ,  when  the 
compressibility  of  matter  differs  up  to  several  orders  of  magnitude.  This 
difference  often  gives  some  unphysical  values  of  variables  appearing  due  to 
errors  of  numerical  approximation  (negative  internal  energy  or  density),  which 
lead  to  a  failing  of  numerical  scheme.  As  it  was  found  in  [16],  some  of 
Godunov-type  second  order  methods  based  on  approximate  Riemann  solver 
are  not  positively  conservative  or  positively  conservative  under  certain 
conditions.  They  can  produce  non-physical  states  with  negative  density  or 
internal  energy  in  rarefaction  waves.  We  based  our  code  on  the  Godunov's 
method  of  numerical  integration  of  Euler  equations  on  a  time  dependent  (so 
called  moving)  grid  [17].  As  it  was  shown  in  [16],  this  scheme  is  positively 
conservative  for  at  least  in  the  case  of  ideal  gas  EOS. 

Let  us  remind  the  main  principles  of  this  code.  In  the  framework  of  this 
approach  only  the  boundaries  of  numerical  grid  are  moved  in  Lagrangian 
fashion.  The  location  of  interior  grid  nodes  are  determined  during  grid 
generation  procedure  by  using  of  coordinates  of  boundary  nodes.  As  the 
material  flows  through  the  grid  cells  there  is  some  diffusion  error  as  well  as  in 
any  Eulerian  method.  Nevertheless  there  exist  several  possibilities  to  reduce 
such  errors  by  generating  "near-Lagrangian"  grids.  In  the  present  version  of 
2D  hydrocode  we  utilize  decomposition  of  numerical  region  onto  subregions 
with  Lagrangian  boundaries,  in  which  grids  are  generated  independently 
using  conformal  mapping  procedure.  This  decomposition  is  provided 
automatically  during  the  computations.  When  every  subregion  consist  of  only 


one  cell.  This  approach  resembles  Lagrangian  method  on  unstructured 
quadrilateral  mesh.  On  the  other  side,  such  a  grid  enables  all  advantages  of  a 
regular  grid.  It  is  known  that  the  numerical  approximation  of  solution  on  an 
orthogonal  grid  is  always  better  than  that  one  on  an  arbitrary  disturbed  grid 
[17].  The  grid  generation  procedure  can  be  governed  by  moving  the 
boundary  nodes  along  the  boundary.  Using  such  a  redistribution  of  the 
boundary  nodes  orthogonal  grids  are  generated  in  subregions.  The  other 
advantage  of  orthogonal  grids  is  a  possibility  to  simplify  a  second  order 
extension  of  Godunov's  scheme. 

In  the  present  report  we  consider  a  second  order  extension  of  moving  grid 
code  and  compare  results  of  calculations  of  hypervelocity  impact  problems 
with  shock-fitting  calculations.  The  strategy  of  shock  fitting  grids  was 
described  in  the  previous  report.  Unfortunately,  the  application  of  this 
approach  is  not  possible  inside  of  numerical  subregions.  It  means,  that  the 
description  of  shock  wave  propagation  and  interaction  inside  of  disturbed 
subregion  can  not  be  done  using  shock  fitting  grids.  That  is  the  reason  why 
the  second  order  extension  of  Godunov's  scheme  is  very  important. 

3. NUMERICAL  ALGORITHM 

In  this  report  we  describe  algorithm  of  computations  only  briefly.  More 
detailed  description  one  can  find,  for  example,  in  [18],  as  well  as  in  the  first 
report.  As  we  have  already  mentioned,  the  numerical  region  is  divided  into 
several  subregions  with  Lagrangian  boundaries  during  the  computations.  Let 
us  consider  the  main  steps  of  computational  algorithm  for  one  subregion. 

The  first  step  of  computations  is  the  displacement  of  the  boundary  of 
subregion.  After  shifting  the  boundary  to  a  new  position  some  segments  of 
the  boundary  can  change  their  type  due  to  interaction  with  boundaries  of  the 
other  subregions.  The  boundary  type  determines  the  boundary  condition 
(what  is  necessary  for  Riemann  problem  initial  data)  and  the  law  of  motion  of 
the  boundary.  For  example,  a  shock-front  boundary  moves  according  to 
Hugein's  principle,  rigid  wall  boundary  does  not  move  and  so  on. 

The  second  step  is  an  iterative  procedure  of  orthogonal  grid  generation 
inside  of  subregion.  This  procedure  is  based  on  the  distribution  of  the 
boundary  nodes.  This  distribution  is  performed  several  times  until  the  grid 
generated  with  conformal  mapping  becomes  orthogonal. 

The  third  step  is  a  solution  of  Riemann  problems  for  inner  zones  and 
calculation  of  fluxes  between  the  neighboring  zones.  One  should  remark  here. 


that  for  an  arbitrary  equation  of  state  the  Riemann  problem  can  be  solved 
only  numerically.  Nevertheless,  we  employ  an  exact  solver  of  this  problem 
only  in  vicinity  of  flow  discontinuity.  The  most  of  Riemann  problem 
computations  are  done  either  in  isentropic  or  in  acoustic  approximation. 

Second  order  extension  of  Godunov's  scheme  can  be  obtained  if  we 
assume  a  piecewise  linear  distribution  of  flow  parameters  inside  of  the  grid 
cells.  To  conserve  the  monotonicity  property  of  Godunov's  scheme  we  use 
the  "minimum  derivative  principle"  proposed  in  [19].  The  main  idea  of  this 
principle  is  to  choose  the  minimum  possible  derivative  when  interpolating  the 
values  from  the  zone  center  to  the  boundary  with  the  neighboring  zone, 
where  the  Riemann  problem  must  be  solved.  If  the  grid  is  orthogonal, 
derivatives  only  in  one  direction  (either  along  the  grid  rows,  or  along  the  grid 
columns)  are  important  when  interpolating  cell-centered  values  to  the 
boundaries  of  the  cells.  This  simplifies  the  realization  of  "minimum  derivative 
principle". 

4.EXAMPLES  OF  CALCULATIONS 

In  the  present  paper  we  consider  an  application  of  our  2D  code  to 
hypervelocity  impact  problem.  One  of  the  most  difficult  task  for  any  2D 
Lagrangian  code  is  to  compute  hypervelocity  penetration  of  projectile  into  a 
thick  target.  We  consider  two  different  approaches  to  numerical  simulation  of 
hypervelocity  impact  of  iron  projectile  moving  with  1 0  km/ s  velocity  on  a 
thick  aluminum  target  . 

First  one  concludes  in  a  usage  of  "near-Eulerian"  grid,  which  covers  the 
region  of  the  target  where  the  shock  wave  propagates. 

The  second  one  is  fitting  of  this  shock  wave  by  moving  the  boundary  of  the 
grid  with  shock  front  velocities  obtained  by  Riemann  solver.  The  second 
method  in  fact,  is  more  accurate  to  resolve  the  shock  front  and  is  more 
economic,  because  it  does  not  require  a  grid  in  undisturbed  regions.  On  the 
other  hand  the  shape  of  numerical  region  is  more  complicated  and,  as  a 
consequence  of  that,  the  grid  generated  is  not  so  orthogonal  as  in  the  first 
case.  The  results  of  computations  corresponding  to  these  cases  together  with 
numerical  grids  are  shown  in  Fig.1.  There  were  two  subregions  used  in  the 
case  of  shock-fitting  grid.  One  of  them  corresponded  to  the  projectile,  the 
other  one  covered  the  disturbed  region  in  the  target.  In  the  other  case  (near- 
Eulerian  grid)  three  subregions  were  used  in  computations. 


Fig.l  Numerical  grid  and  pressure  levels  for  hypervelocity  impact  calculations 
using  shock  fitting  grid  (above  the  axis)  and  "near-Eulerian"  grid  (below  the 
axis). 

It  is  known,  that  first  order  Godunov's  scheme  is  a  very  dissipative  one. 
Our  experience  shows,  that  this  scheme  smears  out  weak  discontinuities  in  the 
case  of  realistic  EOS  much  stronger  than  in  the  case  of  perfect  gas. 
Moreover,  this  dissipation  grows  in  time.  Results  of  computation  of 
hypervelocity  impact  of  iron  projectile  on  aluminum  target  shown  in  fig. 2-3 
confirms  this  fact.  In  contrast  to  the  first  order  scheme,  second  order 
calculations  exibit  nearly  constant  width  of  shock  front  smoothing. 


Fig. 2  Grid  boundaries  and  levels  of  constant  pressure  for  the  first  (below  the 
axis)  and  second  order  accuracy  calculations  of  hypervelocity  impact  of  iron 
projectile  (on  the  right)  on  an  aluminum  target. 


Fig. 3  Grid  boundaries  and  levels  of  constant  pressure  for  the  first  (below  the 
axis)  and  second  order  accuracy  calculations  of  hypervelocity  impact  of  iron 
projectile  (on  the  right)  on  an  aluminum  target. 


As  it  has  been  shown  in  [20]  for  the  case  of  steady  gas  flow,  the  second 
order  extension  of  Godunov's  scheme  is  not  so  sensitive  to  distortions  of 
numerical  grid  as  it's  first  order  counterpart.  We  demonstrate  below  (Fig. 4-5) 
the  results  of  simulation  of  hypervelocity  impact  obtained  using  the  first  and 
the  second  order  Godunov  scheme  both  on  the  shock  fitting  grid  and  on  the 
near-Eulerian  grid.  Our  results  generally  confirm  the  conclusions  obtained  in 
[20].  For  the  second  order  scheme  (Fig. 5)  the  difference  between  the  results 
obtained  on  the  different  grids  is  negligible,  while  for  the  first  order  scheme 
there  is  a  big  difference  in  the  amplitude  of  the  shock  wave  calculated  on  the 
shock  fitting  and  on  near-Eulerian  grid.  It  is  naturally,  that  the  shock  front 
dissipation  is  much  pronounced  in  the  case  of  the  first  order  scheme  than  in 
the  case  of  the  second  order  one. 


Fig. 4  Grid  boundaries  and  pressure  levels  (in  kbars)  for  shock  fitting  grid  and 
"near-Eulerian"  grid.  First  order  scheme. 


Second  order  scheme  tiine=7749  ns 
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Fig. 5  Grid  boundaries  and  pressure  levels  (in  kbars)  for  shock  fitting  grid  and 
near-Eulerian  grid.  Second  order  scheme. 


CONCLUSIONS 

We  have  considered  an  application  of  Godunov's  scheme  on  the  moving 
grids  and  it's  second  order  extension  (Kolgan  scheme)  on  the  example  of 
hypervelocity  impact  calculation.  We  have  shown  that  this  scheme  coupled 
with  the  algorithm  of  decomposition  of  numerical  region  onto  subregions 
shows  a  good  robustness  and  flexibility  for  multimaterial  problems  of 
computational  fluid  dynamics  with  large  distortions.  The  other  advantage  of 
decomposition  of  numerical  region  of  complicated  shape  is  a  possibility  to 
generate  orthogonal  grids  in  subregions 

We  have  developed  second  order  Godunov's  scheme  applying  Kolgan's 
algorithm  of  "minimal  derivative".  It  has  been  found  that  Kolgan's  scheme  is 
not  so  sensitive  to  the  form  of  the  moving  grid  as  original  Godunov's  scheme. 
The  results  of  second  order  calculations  are  much  closer  to  more  accurate 
shock-fitting  calculations. 
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I .  Introduction 


Problems  of  hypervelocity  Impact  attract  attention  of 
scientists  during  the  last  thirty  years  and  occupy  a  significant 
part  In  scientific  literature.  The  urgency  of  Investigations  of 
these  problems  follows  from  the  practice.  It  concerns,  for 
example,  evaluation  of  effectiveness  of  spacecraft  debris 
shields,  analysis  of  asteroid  Impact  events,  problems  of  atomic 
power  stations  safety  and  so  on.  Prom  the  other  side  the 
possibility  of  numerical  simulations  of  these  phenomena  becomes 
real  owing  to  progress  In  theoretical  and  experimental  high 
pressure  high  temperature  physics.  The  success  of  numerical 
simulation  of  hypervelocity  Impact  depends  on  two  factors: 
physical  model  (Including  mechanical,  kinetic, thermodynamic, 
etc.  properties  of  matter)  and  numerical  methods  for  Integration 
of  governing  equations.  The  second  part  of  the  problem  will  be 
discussed  below.  As  to  the  first  one,  the  review  of  various 
substances  models  Is  given,  for  example.  In  [1,  21. 

The  most  widespread  numerical  methods  for  hypervelocity 
Impact  (HVI)  simulation  are  based  on  Langranglan  and  particle 
techniques  [  3-7  ].  Huge  experience  has  been  stored  up  to  now  In 
this  area.  It  becomes  clear,  that  the  main  drawback  of 
Lagranglan  codes  Is  the  restriction  on  deformations  of  flow  and 
that  of  for  the  PIC  codes  Is  the  necessity  of  keeping  of 
’’superfluous”  Information  about  flow  parameters  In  additional 
numerical  grid  cells  Into  which  substance  can  get.  Of  course, 
these  difficulties  are  algorithmic  and  can  be  overcome  by  some 
additional  measures.  For  example,  a  usual  procedure  for 
lagranglan  techniques  Is  re Interpolation  of  flow  parameters 
into  a  new  grid,  when  the  grid  deformations  become  large  [81.  As 
to  Eulerlan  techniques,  one  can  mention  about  tracking  methods 
[9],  that  are  used  to  Increase  the  accuracy  of  calculations  near 
matter-vacuum  and  contact  boundaries  and  to  restrict  the  domain 
of  computations.  Unfortunately  eulerlan  codes  have  one  more 
disadvantage  associated  with  mass  dispersion  problems,  that 
makes  It  difficult  to  calculate  debris  cloud  evolution.  The  are 
several  ways  to  decrease  this  dispersion.  For  example  In  [10] 
both  the  projectile  and  the  shield  were  proposed  to  have  nonzero 
Initial  velocities  moving  towards  each  other  to  minimize 


.advectlon  of  materials  between  cells  In  the  Eulerlan  mesh.  In 
order  to  decrease  dispersion  In  eulerlan  computations  hlgh-order 
accurate  advectlon  schemes  also  were  used  [11,  12]. 

As  both  pure  Lagranglan  and  Eulerlan  codes  have  their  own 
drawbacks,  It  Is  naturally  to  develop  arbitrary  Lagrange-Euler 
(ALE)  codes  combining  the  advantages  of  these  approaches. 
Following  this  Idea,  numerical  simulation  of  penetration  of  a 
hypervelocity  projectile  through  multi-plate  shield  was  done 
[131  using  two  codes:  one  for  computing  of  penetration  stage 
(Eulerlan)  and  another  (Lagranglan)  for  debris  cloud  transport 
calculation.  In  recent  years,  some  algorithms  of  computation 
of  real  matter  dynamics  (Including  elasto-plastlc  and  damaged 
matter  )  developed  previously  only  for  Lagranglan  frame  have 
been  extended  to  the  Eulerlan  frame  of  reference  [12,  14,  151. 
This  opens  new  possibilities  for  HVI  simulations  using  more 
flexible  ALE  codes. 

In  the  present  work  we  used  Godunov  ’  s  scheme  for 
computations  [161.  Godunov’s  method  In  moving  grids  provides 
great  possibilities  and  Includes  both  Lagranglan  and  Eulerlan 
approaches  aS  partial  cases.  It  should  be  marked  that  this 
scheme  has  been  Intensively  explored  for  simulations  of 
aerodynamics  problems,  where  high  accuracy  of  calculations  Is 
necessary.  \ 

In  1968  S.K. Godunov,  A. V. Zabrodin  and  G.P.Prokopov  [17] 
formulated  a  programmatic  approach  to  numerical  Integration  of 
nonsteady  Euler  equations  on  moving  grids.  During  several  years 
the  authors  of  this  report  have  being  working  on  realization  of 
this  approach.  Results  of  this  activity  are  resented  below. 

For  simplicity  the  consideration  Is  given  for  the  Euler 
equations  omitting  of  strength  effects.  Applications  of  this 
approach  to  the  problems  of  high  velocity  Impact  Illustrate 
possibilities  of  2D  code  developed.  As  distinct  from  the  most 
Illustrations  of  possibilities  of  this  method  given  In  [17],  the 
problems  of  hyperveioclty  Impact  are  characterized  by  strong 
deformations  of  matter  and  wide  range  of  matter  state 
parameters:  from  Initially  shock-compressed  or  strongly  heated 
condensed  state  to  finally  expanding  gaseous  or  spalled 
substance.  Large  variations  of  density  lead  to  an  enormous 
Increase  of  geometrical  sizes  of  computed  fields  occupied  by 
the  substance.  Contacts  separating  the  substance  from  vacuum  are 


the  physical  boundaries  of  computational  region.  Therefore  the 
choice  of  Lagranglan  meshes  seems  to  be  a  preferable  strategy 
for  simulations  of  this  class  of  problems.  Unfortunately 
Lagranglan  meshes  crash,  when  the  deformations  become  large, 
because  the  cells  deform,  twist  and  finally  get  tangled.  The 
first  factor  leads  to  great  restrictions  on  calculation  time 
step  and  the  second  makes  further  calculations  Impossible.  As  to 
the  Euler Ian  grid  the  marked  above  uncertainty  In  computational 
region  leads  to  the  necessity  of  keeping  additional  grid  sells 
and  to  corresponding  Increase  of  computational  time.  The  other 
problem  of  Eulerlan  techniques  Is  the  problem  of  contacts 
between  substances. 

Methods  based  on  the  use  of  computational  regions  with 
moving  boundaries  and  generation  Inside  of  these  regions  of 
grids  connected  only  with  the  boundaries  seems  to  unify 
the  advantages  of  Eulerlan  and  Lagranglan  meshes  allowing  to 
overcome  their  drawbacks.  In  the  most  examples  given  below  the 
emphasis  Is  placed  on  demonstration  of  possibilities  of 
application  of  moving  grids  for  treating  the  flows  with  strong 
deformations.  Thus,  for  example,  simulations  of  Impact  problems 
presented  In  this  report  have  been  performed  on  moving  grids 
with  fitting  all  contacts  between  different  media  Including  the 
boundaries  with  vacuum.  During  the  process  of  Interaction  the 
structure  and  configuration  of  boundaries  change:  segments  of 
contacts  can  disappear, for  example  because  of  rebound,  and  new 
segments  can  originate  due  to  Involving  of  new  bodies  Into  the 
process  of  Interaction.  The  flow  region  under  consideration  can 
be  subdivided  Into  several  subregions.  Boundaries  of  subregions 
can  be  both  as  original  flow  singularities  (characteristics, 
boundary  between  substance  and  vacuum,  the  contact  boundary 
between  two  substances,  the  axis  (plane)  of  symmetry,  the  shock 
wave)  and  so  moving  In  accordance  with  their  physical  nature. 
They  can  be  also  Inserted  for  some  computational  reasons. 

As  It  was  mentioned  above,  the  basic  attention  Is  paid  to 
computational  aspects  of  simulation.  Nevertheless,  the  algorithm 
developed  makes  It  possible  using  of  different  physical  models 
In  the  frame  of  Godunov's  scheme.  It  concerns,  for  example,  the 
Rlemann  problem  solver,  designed  for  the  case  of  arbitrary  EOS. 
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II.  Physical  model 


Governing  Equations 

One  fluid  one  temperature  hydrodynamic  model  Is  used  In 
simulations.  Governing  Euler  equations  are  written  In  the  form 
of  conservation  laws  of  mass,  momentum,  energy  and  porous 
volume,  what  Is  necessary  for  Godunov’s  method.  The  conservative 
form  of  the  motion  equations  allows  to  describe  both  as 
continuous  and  so  discontinues  flows. 

^  +  vpu  =0, 

+vpu  ^  u=  -vp(l-Vp), 

^  +  vpuh=  -vpud-Vp)^  (2.1) 

Here  p-  density,  u-  mass  velocity,  h-total  specific  energy, 
p-  pressure,  V^-  porous  volume  per  unit  volume,  (p-  porous 
growth  rate.  For  cf)  we  use  semlemplrlcal  spall  damage  kinetic 
built  on  the  base  of  free  surface  motion  registration  [11. 
Nevertheless  this  model  of  spall  fracture  Is  valid  only  lor 
small  (up  to  several  percents)  porosity.  When  the  porosity  of 
matter  becomes  greater  then  some  value,  the  substance  does  not 
resist  to  expansion.  In  this  case  the  volume  occupied  by  the 
solid  component  does  not  change  anymore  and  the  porosity  grows 
with  the  Increase  of  specific  volume  of  matter. 

Equations  of  state 

Semlemplrlcal  wide-range  equations  of  state  (EOS)  are  used 
[18]  to  determine  the  pressure,  the  temperature  and  the  sound 
velocity  as  functions  of  Internal  energy  per  unit  mass  s  and 
the  density  of  solid  component  p^^  connected  with  density  p  as 
pp=p/(I-Vp).  These  equations  are  built  on  the  base  of 
available  Hugonlot  data  and  provide  for  correct  asymptotes  In 
the  cases  of  extremely  high  energies  and  small  densities. 


III.  Algorithm  of  computations 


1.  Scheme  of  integration  on  the  moving  grid 

Eqs.(2.1)  are  Integrated  on  the  moving  grid  by  Godunov’s 
procedure  [191.  The  scheme  of  the  algorithm  Is  shown  In  Flg.1. 
In  comparison  with  the  traditional  one  ([17])  this  algorithm 
enables  a  possibility  of  the  numerical  region  division  Into  many 
subregions  with  Independent  grid  generation  In  every  subregion. 
An  Important  problem  which  should  be  solved  by  that  Is 
Interaction  between  the  boundaries  belonged  to  different 
subregions . 

2.  Movement  of  boundaries 

The  first  step  of  the  algorithm  Is  the  displacement  of 
subregion’s  boundaries.  The  boundary  of  each  subregion  Is 
represented  as  a  coordinate  array  nodes.  Any  boundary 
between  adjacent  subregions  represents  a  unification  of  nodes 
belonging  to  boundaries  of  these  subregions.  The  coordinates  of 
nodes  are  also  stored  In  special  arrays  establishing  the 
sequence  of  nodes.  Nodes  corresponding  to  the  beginning  and  to 
the  end  of  different  segments  of  boundaries  are  also  marked.  It 
means  that  each  boundary  segment  Is  represented  as  a  sequence  of 
points  belonging  to  the  boundaries  of  different  subregions  and 
that  of  special  end-points  of  boundaries.  The  total  number  of 
nodes  representing  boundary  exceeds,  generally  speaking,  the 
number  of  cells  being  In  contact  from  the  side  of  adjacent 
subregions . 

To  determine  velocity  of  a  boundary  rib,  the  Rlemann 
problem  with  flow  parameters  from  adjacent  cells  being  In 
contact  along  every  link  Is  solved.  In  case  the  Rlemann  solver 
yields  a  negative  pressure  a  rebound  along  the  link  Is  supposed 
to  be  happened.  This  fact  Is  taken  Into  account  by  Inserting 
some  new  special  boundary  points  In  the  array.  In  more  complex 
physical  model  the  coupling  force  can  be  prescribed  to  the 
contact  and  a  rebound  can  be  Inserted  when  the  tense  strengths 
exceeds  the  value  of  this  force. 


Flg.1  Block-scheme  of  the  algorithm 


,  The  velocities  of  nodes  can  he  calculated  after  solving  of 
Rlemann  prohlem  at  the  rib  connecting  the  nodes.  In  the 
Godunov’s  method  the  flow  parameters;  velocity,  density, energy 
are  assumed  to  be  uniform  inside  the  meshes.  For  example,  one 
should  take  flow  parameters  from  the  meshes  A  and  B  to  calculate 
velocity  Uj  and  so  from  the  meshes  A  and  C  to  calculate  velocity 
U^,  see  fig. 2  There  may  be  different  ways  of  nodes  velocities 
Interpolation.  The  linear  Interpolation  proposed  by  Godunov 
[17]  gives  the  formula 


12+  2  1 

u=  -  (  3.1  ) 

1  +1 

1  2 

How  It ' s  seen  from  (  3.1)  the  shorter  the  rib  Is ,  the  more 
contribution  gives  It  In  the  node  velocity.  (For  example.  If 
ii»i2  then  u  s  u^)  From  another  point  of  view.  If  we  want  to 
minimize  the  middle-square  deviation  between  the  boundary 
displacement  obtained  from  the  Rlemann  problem  solution,  when 
every  rib  moves  with  an  appropriate  "big"  velocity,  and  the 
boundary  displacement  determined  by  the  nodes  movement  one 
should  use  the  following  formula  (how  It  has  been  shown  In 
[201) 


*  11+22 
u  =  -  (3.2) 

1  +1 
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In  present  algorithm  we  use  an  average  way  of  Interpolation 
between  (3.1  )  and  (3.2  )  ,  which  leads  to  a  simple  formula 

u=0.5*(u  +u* )=0.5* (u^+u^ )  (3.3  ) 

The  formula  (3.3)  combines  the  advantage  of  (3.3)  with 
simplicity. 


In  this  version  of  moving  grid  code  a  procedure  of  boundary 
node  location  correction  Is  proposed,  which  allows  to  generate 
orthogonal  grids  In  subregions.  For  this  purpose,  boundary  nodes 
are  displaced  along  the  boundary  at  any  time  step  to  ensure 
orthogonality,  between  the  boundary  and  correspondent  Inner  grid 
line.  Nodes  shift  Is  limited  by  settled  maximum  ratio  between 
ad;]acent  the  node  rib's  lengths  In  our  calculations  the 
maximum  ratio  was  1.1.  In  order  to  not  distort  the  shape  of 
numerical  region  during  nodes  distribution  the  node  Is  shifted 
along  smooth  curve  line  which  approximates the  boundary.  The  'z 
curvature  of  the  line  Is  changed  from  node  to  node  using 
line  Interpolation  formula  for  curvature  radius.  After  nodes 
displacement  new  grid  Is  generated  In  subregions. 

Grid  generation  algorithm  Is  based  on  a  univalent  mapping 
of  a  rectangular  region  at  a  parametric  plane  (^,r])  onto  a 
subregion  under  consideration  at  physical  plane  (x,y)  of  flow. 

The  problem  can  be  formulated  mathematically  In  a  following  way. 

The  flow  region  at  plane  (x,y)  Is  considered  to  be  a  curvilinear 
quadrangle  with  tops  A,  B,  C,  D  (fig. 3).  The  univalent  mapping 
of  this  quadrangle  onto  the  rectangular  A%B',C',D'  Is 
determined  by  solution  of  a  system  of  elliptical  equations 


^xx'^'^yy  ^  ’'^xx'^"’^yy  (2*4) 

It  can  be  easily  proved  that  this  mapping  Is  univalent, 
for  example  considering  the  principle  of  maximum  of  harmonic 
function.  In  reality,  to  solve  the  system  of  equations 
relative  to  x,  y  as  independent  variable  parameters  Is  more 
convenient 


a 


a  x^^  -  2b  x^^  +  g  x^  =0  (3.5) 

g  y^  =0 

y^  b  =  +  y^y.^  g  =  +  y| 


using  the  following  boundary  conditions: 


x(e)=x(e,T]i) 

x(r))=x(e(j,'n) 

x(o=x(e,iii.) 

x(ri)=x(eu,T)) 


yco^yce.rii) 

y(T))=y(e(j,r]) 

y(o=y(e,T)i.) 

y(T])=y(eu.T)) 


at  AB  (3.6) 
at  BC 
at  CD 
at  DA 


where 


The  transition  to  discrete  analog  of  boundary  conditions  for 
this  system  of  quaslllnear  elliptical  equations  can  be 


performed  easily.  If  the  the  grid  nodes 

are  given 

^1=1  1=1 , 

. . .  ,n 

(3.7) 

1)^=3  3=1 , 

then  the  boundary  conditions 

. . .  ,m 

are 

(3.8) 

x(l)=x(l,1) 

y(l)=y(l,1 ) 

at  AB 

x(3)=x(1,J) 

y(3)=y(i  ,3) 

at  BC 

x(l)=x(l,m) 

y(l)=y(l,m) 

at  CD 

x(J)=x(n,:l) 

y(3)=y(n,d) 

at  DA 

Different  algorithms  have  been  proposed  In  [17]  for 
solving  this  problem.  Iterations  considered  to  be  the  most 
efficient  algorithm.  We  modified  this  algorithm  treating 
solution  of  corresponding  equations  as  a  result  of  relaxation 
of  solutions  of  nonsteady  equations.  For  this  purpose  ’’non- 
steady”  terms  dx/dt  and  dy/dt  are  Inserted  Into  the  right  sides 
of  equations  for  x  and  y  correspondingly  [211.  The  resulted 
parabolic  system  Is  solved  by  the  method  of  variable  directions 
At  the  first  half  time  step  terms  containing  derivatives  over  h 
are  considered  to  be  given.  At  the  second  half  time  step 
derivatives  over  x  are  taken  from  calculations  at  the  first 
stage.  The  problem  at  the  each  stage  Is  similar  to  the 
one-dlmenslonal  heat  transport  problem.  This  consideration  Is 
similar  to  Iterations,  but  as  any  relaxation  method  It  provides 
more  possibilities  for  the  choice  of  "time-step”  or,  what  Is  the 
same,  of  relaxation  parameter.  In  particular.  It  can  be  changed 
from  one  grid  node  to  another  to  accelerate  the  convergence  of 
the  method.  The  proposed  algorithm  has  demonstrated  Its 
reliability.  Nevertheless  It  falls  for  boundaries  of  a 
complex  configurations.  In  order  to  overcome  problems  associated 
with  a  grid  generation  the  algorithm  foresees  the  possibility  of 
subdividing  the  flow  region  Into  subregions  with  reasonably 
regular  boundaries. 


REGION  2 


RE  0  I  ON  1 


Fig. 2  Boundary  between  the  dllferent  subregions  (solid  line) 
-"big”  velocities  obtained  from  the  Rlemann 
problem  solution  between  the  meshes  A  (region  1 )  and  B 
(region  2)  and  so  A  and  C,  1^,  1^  -the  lengths  ol 
the  ribs 


3.  Bouruiary's  type  checking 

Using  or  moving  grid  makes  It  much  more  easier  setting  of 
boundary  condition  In  comparison  with  Elerlan  grid.  In 
accordance  with  the  boundary  condition  some  type  ol  boundary  has 
been  assigned  to  every  rib  belonged  to  the  boundary.  There  are 
live  boundary  types  that  have  been  considered:  "Iree  boundary" 
(boundary  with  vacuum),  "contact"  (boundary  with  another 
matter),  "axes"  (axes  of  symmetry),  "hard  wall"  and  so  called 
"special  boundary"  through  which  a  matter  can  flow.  It’s  clear, 
that  the  type  of  the  boundary  can  be  changed  from  one  time  step 
to  another.  Therefore  a  procedure  checking  the  boundary  type  at 
each  time  step  has  been  developed. 
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4.  Nocbs  distribution  along  the  boundary  and  grid 
generation  technique 

Experience  of  the  numerical  simulation  on  the  moving  grid 
shows  that  some  procedure  of  the  nodes  distribution  along  the 
boundary  Is  required  to  generate  a  good  grid  Inside  of  the 
region.  Interpolation  ol  the  boundary  plays  an  Important  role 
for  displacement  of  nodes  along  the  boundary.  How  It  has  been 
noticed  [18],  the  sharp  angels  of  the  boundary  are  smoothed 
from  one  time  step  to  another,  when  the  nodes  move  along  the 
straight  lines  connecting  them.  To  prevent  this  smoothing  on  the 
boundary,  we  approximate  It  by  arcs  linking  the  boundary  nodes. 

In  the  case  of  strong  deformations  of  numerical  region  the 
law  of  nodes  distribution  along  the  boundary  Influences  In  a 
great  power  on  the  grid  generated  Inside  the  region.  Previously 
[20]  the  law  of  the  nodes  distribution  along  the  boundary  was 
chosen  anytime  In  accordance  with  the  specifics  of  the  problem 
to  be  solved.  Such  a  way  required  at  least  several  attempts  to 
find  the  appropriate  law  of  nodes  distribution  for  problem 
considered. 


Fig. 3  The  scheme  of  the  grid  mapping 


,  5.  Riemm  problem  solution  and  flow  parameters 

recalculation  for  next  time  level 

The  basic  point  of  numerical  Integration  of  governing 
equations  Is  the  Rlemann  problem,  which  Is  not  self  similar  If 
any  kinetics  phenomena  Is  considered  (for  example,  stress 
relaxation  or  mlcroporous  growth).  In  the  model  used  the  porous 
growth  Is  also  Included  Into  consideration.  Therefore  we  take  It 
to  be  ’’frozen”  when  the  time  of  the  fracture  kinetics  t  Is 

k 

much  more  than  the  time  step  At  and  on  the  contrary.  Hugonlot 
adlabates  and  Isoentrops  for  solid  matter  are  taken  for  solving 
Rlemann  problem  In  the  first  case  and  those  for  porous  matter 
are  used  In  the  second  case.  This  approach  has  been  previously 
used  [22]  In  calculations  of  vlsco-elastlc  flows  by  Godunovs 
method. 


J>G 

contact 


Fig. 4  Discontinuity  decay:  1  and  2  designate  undisturbed 
regions,  1 ’  and  2’  -  regions  behind  the  waves, 
dashed  line  -  contact  discontinuity 


Fig. 5  U-P  plane.  Points  1  and  2  designate  Initial  flow 
parameters  from  the  left  and  from  the  right  side 
of  the  rib.  II  means  the  normal  to  rib  velocity. 

Unlike  the  perfect  gas.  In  the  case  of  an  arbitrary  EOS 
the  Rlemann  problem  can’t  be  solved  exactly  by  Iteration  scheme 
without  numerous  EOS  evaluations.  One  should  mention  here  about 
some  approximate  Rlemann  solvers  [23,241  that  can  be  useiul  also 
for  the  problems  considered.  Initial  flow  discontinuity  decays 
Into  the  contact  discontinuity,  the  left  wave  and  the  right 
wave,  those  may  be  the  shock  wave  or  rarefaction  wave  depending 
on  the  Initial  parameters  (see  Flg.4).  The  mass  velocities  and 
the  pressures  behind  the  right  wave  (2-2 ’ )  are  equal  those 
behind  the  left  wave  (1-1’).  These  so  called  ”blg  parameters” 
(PG  and  UG)  can  be  determined  as  Intersection  of  left  and  right 
wave  velocities  functions  In  U-P  plane  (Fig. 5  ). 

Densities  and  energies  In  the  regions  1’  and  2’  can  be 
calculated  If  velocity  UG  and  pressure  PG  are  known.  Velocity 
functions  FL  and  FR  are  built  using  a  real  equation  of  state.  To 
build  the  simple  wave  velocity  function  (2 ’-2)  sound  velocity  as 
a  function  of  pressure  and  density  Is  sufficient.  In  the  case  of 
the  shock  wave  (1-1’)  the  energy  as  a  function  of  pressure  and 
density  Is  necessary.  In  order  to  calculate  the  fluxes  of  mass, 
energy,  porous  volume  and  momentum  the  rib  velocity  has  to  be 
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■known.  Then,  the  Iluxes  can  be  determined  choosing  the  How 
parameter  from  the  region  where  the  rib  velocity  Is  located 
(region  2’  In  the  case  shown  In  Fig. 4  ).  The  velocity  of  the  rib 
V  can  be  determined  11  the  locations  of  the  rib  at  the 
neighboring  time  levels  are  known  (see  Flg.6).  Generally 
speaking,  the  surface  covered  by  the  rib  during  It’s  movement  Is 
not  the  plane.  For  simplicity  we  replace  It  by  the  plane 
following  the  way  proposed  by  Godunov  [171.  In  accordance  with 
this  method  the  normal  component  of  the  rib  velocity  Is 
calculated  as 

b  b’  a’ 

- T-aTT,  <3. 9) 

-r  t  +  A  t  /  2 
ab 

Here  A  -  Is  the  area  covered  by  the  rib  ab  during  the  time 
step  At  in  the  plane  XY,  L  .-the  length  of  the  rib  at  the  time 
t+At/2. 


Flg.6  Determination  of  rib  velocity 
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6.  Conclusion  remark 

Naturally  the  description  of  the  algorithm  given  above  Is 
at  least  schematic  and  except  of  some  modifications  and  details 
represents  only  the  Interpretation  of  Ideas  [171.  However  we 
considered  the  program  approach  as  being  already  given  and  the 
main  problem  Is  to  put  It  Into  practice  for  specific  problems  of 
physics.  On  this  way  one  can  face  with  a  lot  of  concrete 
problems  of  various  extent  of  complexity:  mathematical, 
physical,  algorithmic.  One  of  the  most  Important  points  Is  an 
experience  of  working  with  given  complex  of  program.  We  are  sure 
that  the  difficulties  connected  with  realization  of  the  approach 
will  be  rewarded  by  Its  huge  possibilities.  The  method  allows  at 
least  In  principle  to  fit  flow  singularities  :  shocks  , 
contacts,  characteristics.  From  this  point  of  view  the  method 
unifies  main  advantages  of  characteristics  methods  with  the 
advantages  of  uniform  computational  schemes.  A  promising  future 
we  see  In  Integration  of  this  method  with  analytical  methods  of 
gas  dynamics.  This  method  provides  a  possibility  to  determine 
dependence  and  Influence  domains  and  to  perform  calculations  at 
this  minimal  domains  area.  From  this  point  of  view  It  can  be 
called  economic. 

Last  years  are  being  marked  by  the  Intense  activity  In 
developing  of  high  order  accuracy  modifications  of  Godunov’s 
methods.  Review  of  these  methods  can  be  found  for  example  In 
[251.  These  modifications  as  a  rule  can  easily  be  Included  Into 
the  body  of  the  original  Godunov’s  algorithm.  Modifications 
concern  the  Rlemann  problem  procedure  and  formulae  for 
recalculation  of  flow  parameters  for  the  new  time  step. 
Development  of  high  order  resolution  schemes  In  moving  grids 
fitting  flow  singularities  Is  of  great  Importance  [261.  They  will 
be  considered  more  detally  In  the  next  report. 
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.IV.  Demonstrative  examples. 

1 .  TEST  1 :  impact  of  a  plane  alminm  plate  onto  a 
semiinfinite  alminm  target 

The  simplest  test  problem  lor  2D  code  Is  an  Impact  ol  a 
plane  pro;Jectlle  on  a  semllnflnlte  target.  In  this  case  the 
problem  becomes  to  be  one-dlmenslonal  and  11  Hugonlot  adlabates 
ol  materials  are  known  the  pressure  and  the  shock  velocity  can 
be  easily  calculated.  The  time  when  rarelactlon  wave  Irom  the 
back  side  ol  the  projectile  reaches  the  shock  Iront  can  be  also 
determined.  Results  ol  a  plane  aluminum-  on-  aluminum  Impact 
with  8  km/s  velocity  are  shown  In  Fig.  7-8.  Region  ol 
computations  consists  ol  two  subregions  corresponding  to  the 
target  and  the  Impactor.  lelt  boundary  ol  the  target  moves 
Initially  with  a  characteristic  velocity  (sound  speed)  Inside 
aluminum  target.  Then  shock  Irom  contact  reaches  this  boundary 
and  It  becomes  moving  as  shock  discontinuity  (see  llg.  8).  This 
transition  Irom  sound  to  shock  regime  ol  boundary  movement  Is 
automatic  because  the  velocity  ol  the  boundary  Is  obtain  Irom 
Rlemann  problem  solution  algorithm.  As  It  can  be  easily  proved, 
numerical  result  Is  agreed  well  with  the  analytical  calculation 
ol  the  shock  pressure  and  velocity. 

2.  TEST  2:  influence  of  the  grid  refinement  on  the 
results  of  nmerical  simulation  of  penetration 
problems 

Calculation  with  dlllerent  levels  ol  accuracy  were 
perlormed  lor  aluminum  sphere  ol  1  cm  diameter  penetrating 
aluminum  plate  ol  2  mm  thickness  at  10.1  km/s.  Numerical  grids 
and  levels  ol  constant  pressure  corresponding  the  time  moment  ol 
710  ns  are  shown  In  llg.9-11.  In  the  case  presented  In  llg.  11 
accuracy  was  two  times  more  than  that  In  the  llg.  10  and  lour 
times  more  than  In  the  llg. 9.  It  Is  seen  that  the  higher  the 
accuracy,  the  better  Is  shock  Iront  resolution  In  the 
projectile.  More  line  grid  allows  also  to  resolve  the  low 
density  vaporized  matter  blowing  oil  the  target  during  Impact, 
At  the  lollowlng  stage  ol  penetration  shown  In  llg.  12-1 3  line 
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-  grid  gives  less  thickness  of  the  target ’ s  substance 
surrounded  debris  of  the  projectile.  And  of  course  as  a 
consequence  of  higher  pressure  there  Is  a  higher  expansion 
velocity  of  debris  cloud  In  this  case. 

3.  Interaction  of  spherical  projectile  with  two  layered 
alminm  target 

Result  of  computation  of  spherical  aluminum  projectile 
Interaction  with  a  spaced  two  layered  aluminum  target  Is 
presented  In  fig. 14.  Impact  velocity  for  the  problem  considered 
Is  8  km/s.  One  should  note  that  not  more  than  1000  meshes  Is 
used  In  this  computation.  The  number  of  meshes  grows  during 
computation  to  ensure  some  fixed  level  of  accuracy.  Criteria  for 
grid  refinement  Is  not  pure  geometric.  When  pressure  gradient 
and  density  of  matter  decrease  some  coalescence  of  meshes  takes 
place. 

Using  moving  grids  demonstrates  a  substantial  economy  of 
reserved  computer  memory  In  comparison  with  a  fixed  eulerlan 
grid.  One  can  calculate  that  about  30  000  -  40  000  meshes  Is 
required  to  perform  computation  of  this  problem  using  eulerlan 
grid  with  the  same  accuracy. 

4.  Hypervelocity  impact  of  iron  projectile  on  a  target 
with  a  conical  hole 

This  example  demonstrates  possibilities  of  moving  grid 
adaptations  to  region’s  form  In  the  case  of  strong  deformations. 
In  the  beginning  of  Interaction  shock  pressure  of  about  0.5  Mbar 
Is  generated  In  PMMA  disk,  while  pressures  In  Iron  target  reach 
about  2  Mbar.  Then  two  Jets  of  liquid  PMMA  are  formed  moving 
along  the  walls  of  conical  channel  (fig.  16).  Collision  of  the 
Jets  on  the  axis  leads  to  cumulation  effect.  Pressure  behind  the 
projectile  reaches  to  2  Mbar  (fig. 17).  After  Impacting  with  the 
Iron  projectile  these  Jets  breach  between  the  projectile  and  the 
walls  of  the  channel  accelerating  projectile  (fig.  18). 
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Fig. 8  Pressure  profiles  at  different  time  moments 


Fig. 9  Grid  and  pressure  levels.  Fig. 10  Grid  and  pressure  levels 
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Fig. 11  Grid  and  pressure  levels  (best  accuracy) 


step  100  kbar,  time:  710  ns 
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Fig. 12  Grid  and  velocity  levels:1-10.6,  2-10.3,  3-9.9,  4-9,6, 
5-9.2,  6-8. 8, 7-8. 5,  8-8.1,  9-7.7,  10-7.4  km/s 
Time:  1900  ns 
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Fig. 13  Grid  and  velocity  levels  (the  same  as  for  rig.i2) 
. . .  .  Time:  1900  ns 
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Fig. 14  Grids  and  density  levels  at  dlflerent  time  moments. 


Fig. 17  Grid  and  pressure  levels 


Fig. 18  Grid  and  pressure  levels 
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Introduction 


One  of  the  most  important  problem  of  numerical  solution  of  nonlinear 
hyperbolic  conservation  equation  is  to  obtain  high  accuracy  of  the  solution  in 
both  discontinous  and  continuos  regions.  In  recent  years,  a  lot  of  so  called 
high  resolution  schemes  have  been  proposed,  for  example,  second  order  TVD 
scheme  of  Roe  [1,2]  and  Sweby  [3];  upwind  TVD  of  Harten  [4]  and  Yee  et 
al.  [5];  symmetric  TVD  of  Chakravarthy  and  Osher  [6];  ENO  scheme  of  Harten 
et  al.  [7,8],  FCT  scheme  of  Zalesak  [9]  McDonald  and  Ambrosiano  [10]; 
MUSCL  schemes  of  van  Leer[11],  Goodman  and  LeVeque  [12],  Davis  [13], 
Colella  [14],  and  PPM  of  Colella  and  Woodward  [15].  The  most  of  schemes 
mentioned  above  were  originally  derived  for  the  case  of  perfect  gas  and  for 
ID  geometry  (or  2D  for  fixed  Eulerian  grids).  For  numerical  simulations  of 
hypervelocity  impact  problems,  when  realistic  physical  models  are  used,  we 
face  with  additional  problems.  This  is,  for  example,  a  variety  of  matter 
properties  in  different  phase  states  (solid,  liquid,  gas,  plasma)  ,  when  the 
compressibility  of  matter  differs  up  to  several  orders  of  magnitude.  This 
difference  often  gives  some  unphysical  values  of  variables  appearing  due  to 
errors  of  numerical  approximation  (negative  internal  energy  or  density),  which 
lead  to  a  failing  of  numerical  scheme.  As  it  was  found  in  [16],  some  of 
Godunov-type  second  order  methods  based  on  approximate  Riemann  solver 
are  not  positively  conservative  or  positively  conservative  under  certain 
conditions.  They  can  produce  non-physical  states  with  negative  density  or 
internal  energy  in  rarefaction  waves.  We  based  our  code  on  the  Godunov's 
method  of  numerical  integration  of  Euler  equations  on  a  time  dependent  (so 
called  moving)  grid  [17].  As  it  was  shown  in  [16],  this  scheme  is  positively 
conservative  for  at  least  in  the  case  of  ideal  gas  EOS. 

Let  us  remind  the  main  principles  of  this  code.  In  the  framework  of  this 
approach  only  the  boundaries  of  numerical  grid  are  moved  in  Lagrangian 
fashion.  The  location  of  interior  grid  nodes  are  determined  during  grid 
generation  procedure  by  using  of  coordinates  of  boundary  nodes.  As  the 
material  flows  through  the  grid  cells  there  is  some  diffusion  error  as  well  as  in 
any  Eulerian  method.  Nevertheless  there  exist  several  possibilities  to  reduce 
such  errors  by  generating  "near-Lagrangian"  grids.  In  the  present  version  of 
2D  hydrocode  we  utilize  decomposition  of  numerical  region  onto  subregions 
with  Lagrangian  boundaries,  in  which  grids  are  generated  independently 
using  conformal  mapping  procedure.  This  decomposition  is  provided 
automatically  during  the  computations.  When  every  subregion  consist  of  only 


one  cell.  This  approach  resembles  Lagrangian  method  on  unstructured 
quadrilateral  mesh.  On  the  other  side,  such  a  grid  enables  all  advantages  of  a 
regular  grid.  It  is  known  that  the  numerical  approximation  of  solution  on  an 
orthogonal  grid  is  always  better  than  that  one  on  an  arbitrary  disturbed  grid 
[17].  The  grid  generation  procedure  can  be  governed  by  moving  the 
boundary  nodes  along  the  boundary.  Using  such  a  redistribution  of  the 
boundary  nodes  orthogonal  grids  are  generated  in  subregions.  The  other 
advantage  of  orthogonal  grids  is  a  possibility  to  simplify  a  second  order 
extension  of  Godunov's  scheme. 

In  the  present  report  we  consider  a  second  order  extension  of  moving  grid 
code  and  compare  results  of  calculations  of  hypervelocity  impact  problems 
with  shock-fitting  calculations.  The  strategy  of  shock  fitting  grids  was 
described  in  the  previous  report.  Unfortunately,  the  application  of  this 
approach  is  not  possible  inside  of  numerical  subregions.  It  means,  that  the 
description  of  shock  wave  propagation  and  interaction  inside  of  disturbed 
subregion  can  not  be  done  using  shock  fitting  grids.  That  is  the  reason  why 
the  second  order  extension  of  Godunov's  scheme  is  very  important. 

3.NUMERICAL  ALGORITHM 

In  this  report  we  describe  algorithm  of  computations  only  briefly.  More 
detailed  description  one  can  find,  for  example,  in  [18],  as  well  as  in  the  first 
report.  As  we  have  already  mentioned,  the  numerical  region  is  divided  into 
several  subregions  with  Lagrangian  boundaries  during  the  computations.  Let 
us  consider  the  main  steps  of  computational  algorithm  for  one  subregion. 

The  first  step  of  computations  is  the  displacement  of  the  boundary  of 
subregion.  After  shifting  the  boundary  to  a  new  position  some  segments  of 
the  boundary  can  change  their  type  due  to  interaction  with  boundaries  of  the 
other  subregions.  The  boundary  type  determines  the  boundary  condition 
(what  is  necessary  for  Riemann  problem  initial  data)  and  the  law  of  motion  of 
the  boundary.  For  example,  a  shock-front  boundary  moves  according  to 
Hugein's  principle,  rigid  wall  boundary  does  not  move  and  so  on. 

The  second  step  is  an  iterative  procedure  of  orthogonal  grid  generation 
inside  of  subregion.  This  procedure  is  based  on  the  distribution  of  the 
boundary  nodes.  This  distribution  is  performed  several  times  until  the  grid 
generated  with  conformal  mapping  becomes  orthogonal. 

The  third  step  is  a  solution  of  Riemann  problems  for  inner  zones  and 
calculation  of  fluxes  between  the  neighboring  zones.  One  should  remark  here. 


that  for  an  arbitrary  equation  of  state  the  Riemann  problem  can  be  solved 
only  numerically.  Nevertheless,  we  employ  an  exact  solver  of  this  problem 
only  in  vicinity  of  flow  discontinuity.  The  most  of  Riemann  problem 
computations  are  done  either  in  isentropic  or  in  acoustic  approximation. 

Second  order  extension  of  Godunov's  scheme  can  be  obtained  if  we 
assume  a  piecewise  linear  distribution  of  flow  parameters  inside  of  the  grid 
cells.  To  conserve  the  monotonicity  property  of  Godunov's  scheme  we  use 
the  "minimum  derivative  principle"  proposed  in  [19].  The  main  idea  of  this 
principle  is  to  choose  the  minimum  possible  derivative  when  interpolating  the 
values  from  the  zone  center  to  the  boundary  with  the  neighboring  zone, 
where  the  Riemann  problem  must  be  solved.  If  the  grid  is  orthogonal, 
derivatives  only  in  one  direction  (either  along  the  grid  rows,  or  along  the  grid 
columns)  are  important  when  interpolating  cell-centered  values  to  the 
boundaries  of  the  cells.  This  simplifies  the  realization  of  "minimum  derivative 
principle". 

4.EXAMPLES  OF  CALCULATIONS 

In  the  present  paper  we  consider  an  application  of  our  2D  code  to 
hypervelocity  impact  problem.  One  of  the  most  difficult  task  for  any  2D 
Lagrangian  code  is  to  compute  hypervelocity  penetration  of  projectile  into  a 
thick  target.  We  consider  two  different  approaches  to  numerical  simulation  of 
hypervelocity  impact  of  iron  projectile  moving  with  10  km/s  velocity  on  a 
thick  aluminum  target  . 

First  one  concludes  in  a  usage  of  "near-Eulerian"  grid,  which  covers  the 
region  of  the  target  where  the  shock  wave  propagates. 

The  second  one  is  fitting  of  this  shock  wave  by  moving  the  boundary  of  the 
grid  with  shock  front  velocities  obtained  by  Riemann  solver.  The  second 
method  in  fact,  is  more  accurate  to  resolve  the  shock  front  and  is  more 
economic,  because  it  does  not  require  a  grid  in  undisturbed  regions.  On  the 
other  hand  the  shape  of  numerical  region  is  more  complicated  and,  as  a 
consequence  of  that,  the  grid  generated  is  not  so  orthogonal  as  in  the  first 
case.  The  results  of  computations  corresponding  to  these  cases  together  with 
numerical  grids  are  shown  in  Fig.1.  There  were  two  subregions  used  in  the 
case  of  shock-fitting  grid.  One  of  them  corresponded  to  the  projectile,  the 
other  one  covered  the  disturbed  region  in  the  target.  In  the  other  case  (near- 
Eulerian  grid)  three  subregions  were  used  in  computations. 
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Fig.1  Numerical  grid  and  pressure  levels  for  hypervelocity  impact  calculations 
using  shock  fitting  grid  (above  the  axis)  and  "near-Eulerian"  grid  (below  the 
axis). 

It  is  known,  that  first  order  Godunov's  scheme  is  a  very  dissipative  one. 
Our  experience  shows,  that  this  scheme  smears  out  weak  discontinuities  in  the 
case  of  realistic  EOS  much  stronger  than  in  the  case  of  perfect  gas. 
Moreover,  this  dissipation  grows  in  time.  Results  of  computation  of 
hypervelocity  impact  of  iron  projectile  on  aluminum  target  shown  in  fig. 2-3 
confirms  this  fact.  In  contrast  to  the  first  order  scheme,  second  order 
calculations  exibit  nearly  constant  width  of  shock  front  smoothing. 


Fig. 2  Grid  boundaries  and  levels  of  constant  pressure  for  the  first  (below  the 
axis)  and  second  order  accuracy  calculations  of  hypervelocity  impact  of  iron 
projectile  (on  the  right)  on  an  aluminum  target. 


Fig. 3  Grid  boundaries  and  levels  of  constant  pressure  for  the  first  (below  the 
axis)  and  second  order  accuracy  calculations  of  hypervelocity  impact  of  iron 
projectile  (on  the  right)  on  an  aluminum  target. 


As  it  has  been  shown  in  [20]  for  the  case  of  steady  gas  flow,  the  second 
order  extension  of  Godunov's  scheme  is  not  so  sensitive  to  distortions  of 
numerical  grid  as  it's  first  order  counterpart.  We  demonstrate  below  (Fig. 4-5) 
the  results  of  simulation  of  hypervelocity  impact  obtained  using  the  first  and 
the  second  order  Godunov  scheme  both  on  the  shock  fitting  grid  and  on  the 
near-Eulerian  grid.  Our  results  generally  confirm  the  conclusions  obtained  in 
[20].  For  the  second  order  scheme  (Fig. 5)  the  difference  between  the  results 
obtained  on  the  different  grids  is  negligible,  while  for  the  first  order  scheme 
there  is  a  big  difference  in  the  amplitude  of  the  shock  wave  calculated  on  the 
shock  fitting  and  on  near-Eulerian  grid.  It  is  naturally,  that  the  shock  front 
dissipation  is  much  pronounced  in  the  case  of  the  first  order  scheme  than  in 
the  case  of  the  second  order  one. 


Fig. 4  Grid  boundaries  and  pressure  levels  (in  kbars)  for  shock  fitting  grid  and 
"near-Eulerian"  grid.  First  order  scheme. 


Fig. 5  Grid  boundaries  and  pressure  levels  (in  kbars)  for  shock  fitting  grid  and 
near-Eulerian  grid.  Second  order  scheme. 


CONCLUSIONS 

We  have  considered  an  application  of  Godunov's  scheme  on  the  moving 
grids  and  it's  second  order  extension  (Kolgan  scheme)  on  the  example  of 
hypervelocity  impact  calculation.  We  have  shown  that  this  scheme  coupled 
with  the  algorithm  of  decomposition  of  numerical  region  onto  subregions 
shows  a  good  robustness  and  flexibility  for  multimaterial  problems  of 
computational  fluid  dynamics  with  large  distortions.  The  other  advantage  of 
decomposition  of  numerical  region  of  complicated  shape  is  a  possibility  to 
generate  orthogonal  grids  in  subregions 

We  have  developed  second  order  Godunov's  scheme  applying  Kolgan's 
algorithm  of  "minimal  derivative".  It  has  been  found  that  Kolgan's  scheme  is 
not  so  sensitive  to  the  form  of  the  moving  grid  as  original  Godunov's  scheme. 
The  results  of  second  order  calculations  are  much  closer  to  more  accurate 
shock-fitting  calculations. 
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i  Introduction  ; 

In  this  report  we  present  applications  of  previously  described  moving  grid  Godunov  code  to 
several  hypervelocity  impact  problems.  This  method  employs  the  exact  solution  of  Riemann 
problem,  therefore  it  is  well  suited  for  construction  of  "nonhomogeneous"  computational 
algorithms  fitting  flow  singularities.  It  is  clear,  that  to  fit  all  singularities,  which  can  arise  in  the 
flow  is  rather  difficult  algorithmic  problem  esspecially  in  multidimentional  case.  That  is  why  the 
most  of  widespread  computational  methods  of  shock  wave  dynamics  are  based  on 
homogeneous  algorifiUns  treating  in  the  same  fashion  both  continuos  and  discontinuos  domains 
I '  of  flow.  :  '  - 

To  better  resolve  discontinuities  so  called  high  resolution  shock-capturing  schems  are  used. 
Most  of  them  are  second  order  modifications  of  Godvmov  scheme  [1],  which  use  approximate 
Riemann  problem  solvers  [2-5].  One  should  mention  in  this  context  so  called  ENO  schemes 
(essentially  non-oscillatory  schemes)  which  allows  one  to  increase  the  order  of  accuracy  up  to 
fourth  or  even  more  [6-7].  Of  pourse,  the  higher  is  the  order  of  accuracy,  the  more 
complicated  is  the  computational  algorithm.  For  example,  as  it  was  found  in  [8],  the  cost  of  the 
third-order  ENO  scheme  is  four  times  that  of  the  second-order  one.  It  may  appear  the  second 
order  schemes  to  be  more  profitable  when  using  refined  mesh  than  the  fourth-order  ones. 

We  have  chosen  another  approach  to  construction  of  high  resolution  shock-wave  codes, 
namely,  Godunov  scheme  on  a  moving  grid.  The  main  principles  of  this  approach  were 
described  in  [9-10]  and  in  the  previous  HEDRC  report.  Note  that  in  the  framework  of  moving 
grid  code  shock  waves,  contacts  between  different  materials  and  matter-vacuum  boundaries, 
can  be  easily  fitted  .  By  that  the  grid  boundaries  coinsided  with  these  discountinuities  are 
moved  in  accordance  with  the  solution  of  Riemann  problem. 

In  the  present  report  we  demonstrate  the  advantages  of  moving  grid  approach  for  solving 
four  different  problems. 

The  first  one  is  an  impact  of  thin  plane  projectile  on  a  thin  screen.  A  great  difference  of 
longituidal  and  transverse  dimentions  leads  to  a  great  number  of  zones  when  using  a 
Lagrangian  grid,  which  covers  all  the  projectile  and  the  target.  Our  method  allows  one  to 


compute  only  those  domains  where  the  flow  is  two-dimentional,  that  gives  a  substantioi 
economy  of  computer's  resources. 

The  second  example  is  the  numerical  simulation  of  penetration  of  a  cosmic  body  moving  w 
hypervelocity  into  gaseous  atmosphere  with  exponentionally  increasing  density.The  typic 
feature  of  this  problem  is  the  presence  of  various  space  scales.  For  examples,  the  size  of  tl 
body  is  much  larger  than  the  thickness  of  compressed  gas  layer  at  the  face  side  of  the  body  ai 
on  the  other  side  it  is  much  less  than  the  distance  which  this  body  passes  before  it  explodes. 

The  third  problem  is  a  debris  cloud  simulation  produced  by  ball-plate  impacts.  In  th 
problem  the  ball  and  the  plate  have  thicknesses  of  the  same  order,  that  makes  it  difficult  to 
employ  any  Lagrangian  code  because  of  the  great  deformations  of  material.  On  the  other  sid 
it  is  known  that  Eulerian  codes  give  a  great  dispersion  of  mass  due  to  averaging  of  flow  valu 
over  the  fixed  zone  intervals.  Application  of  the  moving  grig  code  in  which  the  computatior 
grid  moves  with  material  near  the  boudaries  and  has  a  minimal  advection  error  in  the  ceni 
grid  region  allows  one  to  overcome  these  diffuculties. 

The  last  example  is  a  numerical  simulation  of  detonation  of  charges  of  finite  diameters.  It 
well  known,  that  a  steady  detonation  can  exist  only  if  the  size  of  the  explosive  is  greater  th; 
some  minimum  value.  To  find  this  size  for  a  charge  of  some  specified  geometry  numeroi 
numerical  computations  are  required.  When  this  size  is  close  to  the  critical  one  the  time  whi< 
is  nessesary  for  setting  up  a  quasi  steady  detonation  regime  is  very  long.  It  means  the  we  mu 
be  able  to  track  the  propagation  of  the  front  over  a  long  distance.  As  the  region  i 
computation  increases  with  time,  a  grid  refinement  is  required  to  maintain  some  specified  lev 
of  computation  accuracy. 

Numerical  simulation  of  impact  of  projectiles  with  a  great  D/L  ratio 

Impactors  with  a  great  D/L  ration  (where  D  is  the  diameter  and  L  is  the  thickness)  are  usee 
in  shock-wave  experiments  in  order  to  simplify  their  interpretation  reducing  the  last  one  to  11 
case.  Nevertheless,  the  flow  becomes  two  dimentional  when  the  rarefaction  waves 
from  the  periphery  reach  the  center  of  the  projectile.  To  account  2D  effects  at  the  late  stages 
interaction  2D  simulation  is  nessesary.  In  the  present  report  we  demonstrate  one  example  of 


application  of  moving  grid  to  effective  solution  of  such  a  problem.  In  fig.l(a-d)  numerical  gri 
and  levels  of  constant  density  are  presented  for  high  velocity  impact  of  1mm  thick  aluminum 
disk  on  2mm  thick  aluminum  foil  at  8  km/s  velocity.  The  diameter  of  the  disk  is  4  cm 
In  the  beginning  of  computation  we  specify  two  small  numerical  regions.  One  of  them 
represents  the  edge  of  the  disk  and  the  other  corresponds  to  the  fragment  of  the  foil,  which  i^ 
in  contact  with  this  edge.  Then  we  mov#  the  boundary  segment  separating  the  ID  and  2D  flo\ 
regions  with  characteristics  velocity  towards  to  the  axis  of  simmetry.  The  other  fragments  of 
the  boundary  are:  free  boundaries  (back  sides  of  the  disk  and  the  foil),  contact  between  the 
disk  and  the  foil  and  the  front  of  the  shock  wave  runing  away  from  the  center.  Rarefaction 
waves  coming  from  the  back  surfaces  of  the  projectile  and  the  target  lead  to  an  expantion  of 
shocked  material  and  to  generation  of  tensile  stresses.  That  is  why  the  contact  between  tlie  di 
and  the  target  disappears  .  The  type  of  the  boundary  segments  is  changed  automatically  from 
contact  boundary  to  free  one  (fig.  lb)  For  the  later  times  the  size  of  the  flow  under 
consideration  is  much  more  than  initial  one  (see  fig.  Id).  The  straight  line  fragment  ot  the 
boundary  separates  the  flow  region  where  the  transversal  rarefaction  waves  are  important. 

In  addition  to  this  example  one  should  note,  that  there  is  a  lot  of  problems  containing 
different  time  and  space  scales  which  could  be  effectively  solved  using  moving  grid  code,  for 
example,  investigation  of  thin  foil  acceleration  with  laser  and  ion  beams,  perforation  of 
multilayered  spaced  shields  and  so  on. 

Computation  of  large  asteroid  penetration  in  Jupiter's  atmosphere 

The  problems  of  asteroid  hazard  are  of  great  interest,  because  the  number  of  ateroids  movii 

,•  if  . 

in  vicinity  of  Earth  is  increasing  with  time.  An  impact  of  pereodic  comet  Shoemaker-Lev) 
on  Jupiter  took  place  recently  (in  July  1994).  The  consequences  of  this  impact  were  observt  : 
both  from  the  Earth  and  space  sattelits.  In  this  report  we  demonstrate  an  example  of  numeric;!! 
simulation  of  the  penetration  of  one  comet  fragment  of  1  km  diameter  into  Jupiter's 
atmosphere.  The  asteroid  was  approximated  by  a  spherical  incompressible  body  moving  with  a 
velocity  of  60  km/s.  Due  to  its  high  velocity  it  passses  a  distance  much  more  than  its  diametei 
before  it  explodes.  As  the  density  of  the  atmosphere  grows  with  distance  the  flow  is  unsteady. 


Extremly  high  velocity  of  impact  leads  to  heating  and  ionization  of  a  gas  flow  after  shock 
compression.  Since  the  flow  behind  the  shock  fi-ont  is  responsible  for  ionization  and  luminosity 
of  gas,  which  can  be  observed  an  accurate  resolution  in  this  region  is  required.  It  is  clear,  that 
application  of  Eulerian  codes  to  this  problem  will  lead  to  a  great  difiiision  of  shock  front  and 
contact  interfaces.  Eulerian  computations  also  require  a  great  number  of  computational  zones 
to  cover  the  regions  in  which  the  asteroid  moves.  That  is  why  Eulerian  computations  of  this 
problem  were  performed  in  'reverse  ballistic'  sense  using  an  atmosphere  moving  towards  an 
initially  stationary  fragment  [11]. 

To  avoid  computations  in  undisturbed  domains  we  surrounded  the  asteroid  region  by  a  tli 
region  of  the  gas,  whose  boundary  moves  with  the  characteristics  velocity.  The  size  of  the  e 
region  grows  in  time  but  we  can  exclude  from  consideration  those  parts  of  the  gas  regioK, 
which  are  far  from  the  asteroid  and  are  not  of  our  interest.  The  results  of  computation  are 
shown  in  fig.2-4.  In  the  beginning  (fig.2)  of  penetration  the  gas  detaches  from  the  back  side  of 
asteroid  and  a  gas  jet  appears  moving  in  the  opposide  direction.  Then  we  cut  off  the  part  of  gri  .s 
region  corresponding  to  this  jet  (fig.3).  Fig.4  presents  the  results  of  simulation  for  the  time 
when  the  asteroid  has  passed  about  80  km  in  atmosphere.  The  pressures  at  the  frontal  surface 
of  the  asteroid  become  about  0.3-0.5  kbar,  that  are  close  to  the  strength  of  material.  Furthei 
computations  require  to  account  for  the  fracture  of  matterial.  In  accordance  with  [121 
fragmentation  of  the  asteroid  leads  to  a  growth  of  the  heat  flux  inside  it,  because  of  signifficai  .! 
increase  of  the  effective  surface  area.  This  gives  rise  to  the  rapid  transformation  of  the  asteroi ' 
material  from  the  condensed  to  gaseous  state.  If  we  suppose,  that  the  fracture  ai 
fragmentation  of  asteroid  is  due  to  the  stress  gradient  at  its  surface,  we  obtain  the  distance  < 
about  100  km,  where  the  explosion  takes  place.  This  agrees  with  recent  observations. 


Numerical  simulation  of  propagation  of  debris  cloud  produced  by 
ball-plate  impact. 

The  results  of  a  computation  of  a  20  g  spherical  lead  projectile  striking  a  lead  plate  is 


presented  in  Fig.5(a,b).  The  impact  velocity  for  this  problem  is  6.6  km/s.  Fig.6  represents  tlie 
experimental  X-ray  photograph  of  the  same  problem  at  the  time  moment  30  p.  m  presented  in 
[13],  Simulation  results  are  in  accordance  with  the  experiment.  The  number  of  zones  used 
grows  during  the  computation  to  ensure  some  specified  level  of  accuracy.  The  criteria  for  grid 
refinement  are  not  purely  geometric.  When  the  pressure  gradient  and  the  density  decrea.^e 
some  coalescence  of  zones  takes  place. 

Using  moving  grids  demonstrates  a  substantial  economy  of  required  computer  memory  in 
comparison  with  a  fixed  Eulerian  grid.  One  can  estimate  that  about  10  000  -  20  000  mesh 
is  required  to  perform  computation  of  this  problem  using  Eulerian  grid  with  the  sai; 
accuracy. 

Simulation  of  detonation  in  liigb  explosive  charges  of  finite  diameter. 

t' 

The  failtu’e  detonation  problem  is  the  problem  of  minimal  explosive  charge  diameter  when 
a  self-sustained  detonation  can  eTdst.  We  present  below  an  example  of  determination  of  this 
diameter  using  numerical  simulation.  The  flow  is  described  by  Euler  equations.  The  only 
difference  with  hypervelocity  impact  problem  is  the  appearance  in  the  energy  equation  of  a 
source  term  which  is  responsible  for  heat  release.  To  check  whether  the  detonation  becomes 
steady  or  not  it  is  necessary  to  calculate  the  evolution  of  detonation  process  along  the  distanc 
of  at  least  several  (perhaps  tens)  diameters.  As  the  downstream  flow  behind  the  shock  is 
determined  by  the  chemical  reaction  kinetics  and,  in  particular,  by  the  value  of  energy 
release,  the  shock  must  be  calculated  accurately. 

The  following  computational  strategy  is  chosen  in  accordance  with  the  problems  listed  above. 
The  computational  region  represents  a  curvelinear  quadrangle.  The  sides  of  the  quadrangle  are: 
the  leading  shock,  the  segment  of  contact  between  detonation  products  and  vacuum,  the 
segment  of  axis  of  symmetry,  and  the  forth  one  is  a  segment  of  straight  line  between  the  free 
surface  and  the  axis  of  symmetry,  which  is  perpendicular  to  the  axis  of  symmetry.  The  velocity 
of  this  side  is  assumed  to  be  directed  along  the  axis  and  to  be  equal  to  max  (u+a,D),  where  i ) 
is  the  shock  velocity  at  the  axis  of  symmetry  and  u+a  is  the  velocity  of  characteristic 


surface  moving  along  the  axis  and  calculated  using  parameters  ofthe  cells  adjacent  to  this 
side.  It  means  that  the  boundary  moves  relative  to  the  matter  ahead  of  it  at  least  with  the  sonic 
velocity.  Therefore  the  flow  parameters,  which  determine  the  fluxes  of  mass,  momentum 
and  energy  throughout  this  side  are  taken  from  internal  cells.  A  pressure  of  200  kbar,  normnl 
density  and  zero  velocity  were  taken  as  an  initial  d^ta  for  calculations.  It  was  found  from 
simulations,  that  the  critical  diameter  is  somewhere  between  2.4  and  3  mm.  Results  of 
simulations  corresponding  to  these  two  cases  are  shown  in  Fig.7(a,b)  In  the  case  of  2.4  mm 
diameter  detonation  decays.  Calculations  for  this  diameter  were  performed  with  different 
initial  pressures.  In  all  the  cases  a  steady  detonation  was  not  obtained.  This  fact  allows  to 
draw  a  conclusion  that  the  failure  diameter  of  TNT  lies  within  the  interval  2.4  and  3  mm. 
This  also  agrees  with  the  experiment  [14]. 

CONCLUSIONS 

We  have  demonstrated  the  robotness  and  flexibility  of  developed  2D  code 
for  hypervelocity  problem  computations.  The  main  advantage  of  this  code 
in  comparison  with  Eulerian  high-order  accuracy  hydrocodes  is  an  accurat< 
treatment  of  multimaterial  interfaces.  On  the  other  side,  arbitrary 
Lagrangian-Eulerian  methods  (ALE),  which  utilize  lagrangian  motion  of 
interfaces  and  permit  an  arbitrary  mesh  motion  inside  the  computation 
region  are  not  always  adapted  to  dynamically  evolving  interface  shape.  It  is 
the  main  reason  to  develope  ALE  codes  on  unstructured  grids  [15]. 

We  overcome  these  difficulties  with  the  help  of  decomposition  of  numerical 
region  onto  subregions.  This  decomposition  is  provided  automatically  during 
the  computation.  By  that  one  can  govern  this  process  excluding  from 
computations  the  regions  which  are  not  of  our  interest. 
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Fig.2.  isjumerical  grid  and  pressure  levels.  Pressure  level  step  is  0.001  kbar, 
maximum  pressure  is  0.01  kbar.  Time  is  55  fis. 


Fig.B.I  Numerical  grid  and  pressure  levels.  Pressure  level  step  is  0.001  kbar, 
maximum  pressure  is  0.014  kbar.  Time  is  80  fis.  Initial  conflguration  is  shown  an 
the  right  side. 
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Fig.4.  Numerical  grid  and  pressure  levels.  Pressure  level  step  is  0.06  kL'ar, 
maximum  pressure  is  0.6  kbar.  Time  is  1.3  s. 
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Brief  description  of  the  code 
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The  language  of  the  code  is  standard  F77  FORTRAN.  So  it  can  be  easily 
installed  on  any  computer  with  f77  compiler.  The  code  consists  of  the 
following  program  modules: 

Function: 

input  of  initial  data  and  transformation 
of  these  data  into  file  "NAM".TMP  for 
continuation  of  calculations  by  module  BIG2T 

carries  out  calculations,  renews  data  for 
continuation  of  calculations  (file  "NAM".TMP) 
and  records  results  of  the  calculation 

rnodification  of  the  task  (possibihty  to 
add  or  delete  subregions  ).  The  changed  data 
are  rewritten  into  the  file  "NAM".TMP 

transformation  of  "NAM".TMP  from  binary  file  to 
ASCII  file  and  back  (it  may  be  useful  for 
transmission  of  calculation  results  from  one 
computer  to  other  one  with  different  binary 
data  presentation) 


Module: 

STARTT- 

BIG2T- 

REP- 

UNFTOF- 


Installation  of  the  code 
To  create  object  modules 

1)  copy  all  files  *.FOR  and  *.INC  to 
one  directory  on  a  hard  disk 

REMARK:  all  *.INC  files  must  be  with  capital  characters 

2)  compile  all  *.FOR  files 

*  before  compilation  probably  may  be  necessary: 

A)  to  set  function  which  determines  CPU  time  in  seconds  in 
SUBROUTINE  ETIME  in  the  file  PC.FOR.  (to  decomment 
correspondent  string  may  be  sufficient) 

B)  to  change  parameters  in  the  file  BIGPAR.INC 
(dimensions  of  arrays) 


3)  create  loading  module  STARTT,  using  object  files: 


startt.o,  blocka.o,  adi.o,  arct.o,  big.o,  contw.o,  univer.o,  datO.o, 
dattst.o,  disto.o,  jup.o,  pc.o,  shr.o,  spect.o,  urst.o,  coef.o 

4)  create  loading  module  BIG2T,  using  object  files: 

basica.o,  blocka.o,  mainxt.o,  spect.o,  tnew2a.o,  urst.o,  adi.o,  arct.o, 
big.o,  buil.o,  builO.o,  buill.o,  build.o,  contw.o,  datO.o,  del.o,  disto.o, 
god2.o,  gmew.o,  jup.o,  pc.o,  ref.o,  shr.o,  split.o,  splitl.o,  univer.o, 
coef.o 


5)  create  loading  module  REP,  using  object  files: 

rep.o,  blocka.o,  spect.o,  urst.o,  adi.o,  arct.o,  big.o,  contw.o,  datO.o, 
del.o,  jup.o,  pc.o,  ref.o,  shr.o,  spHt.o,  splitl.o,  umver.o 

*Note  that  modules  mentioned  above  (STARTT,  BIG2T,  REP)  can  be 
created  using  makefile.  To  generate  first  order  accuracy  code  (BIGIT), 
one  should  replace  tnew2a.*  by  tnewla.*  in  the  makefile. 

Generation  of  a  new  initial  data  Hie 

1)  Chose  one  of  the  existing  initial  data  files  from  the 
file  Ust  below 
TST2.DAT 
VRB1.DAT 
VRB2.DAT 
VRB3.DAT 


2)  Copy  the  chosen  file  into  file  "NAM".DAT,  where  NAM  is 
any  name  consisting  of  1-4  symbols 

(the  names  of  all  the  files  created  during  the  computation 
of  this  problem  will  include  "NAM"  as  a  root  word) 

3)  Correct  contents  of  the  file  "NAM".DAT  with  text  editor 
using  description  of  initial  data  file  presented 

below 


Description  of  the  structure  of  initial  data  file  NAM.dat 

The  initial  data  file  consists  of  two  parts: 

First  part  contains  control  information.This  part  has  the  same  format  for  all 
files  of  this  type.  Let  us  consider  an  example  shown  below,  where  the 
fist  part  of  initial  data  file  TST.DAT  is  presented. 


AL/AL  U=  10.1  KM/C 

JSPR-LOS-LST-IRD-ISP-IRA-ICM-AIR 
0  1  1  2  1  00  1  0 
NBL-MI-ID1-ID2-IK-IGR-ILENPARIVSP 
9900  1  10  10  300  10  -10  6  0 

- CMM— : - H3 - : - HDEL— : - T - STABL— : 

-l.OOOOOE-03  0.00075E-00  O.OlOOOE+00  l.lOOOOE+05  0.50000E+00 


The  first  string  is  comment.  This  is  usually  short 
information  about  the  problem.  (Here  this  is  10.1  km/s  impact 
of  aluminum  on  aluminum).  The  next  strings  consist  of 
title_strings  and  contents_strings.  The  title  strings  show  the 
names  of  variables.  The  contents  strings  are  values  of  these 
variables.  The  assignment  of  the  variables  is  shown  in  the  table 

TABLE 


variables:  meaning: 

ISPR  controls  rebound  of  different  boundary  segments 

0-  calculations  without  rebound 

1-  rebound  for  "contact"  type  boundaries  only 

2-  rebound  for  "contacts  and  rigid  walls"  type 
boundaries 

3-  rebound  for  aU  boimdary  t^^jes 

LOS  1-  cyllindrical  symmetry 

0-  plane  symmetry 

LST  1-  there  are  rigid  walls 

0-  rigid  walls  are  absent 

IRD  0-  to  ignore  isp 

1-  to  continue  calculations  from  the  restart  data 

2-  to  correct  in  dialog  regime  the  number  of  time  steps 
after  which  the  spUtting  of  numerical  regions  into 
subregions  will  be  switched  on 

ISP  1-  to  change  in  dialog  regime  the  type  of "  contact " 
between  subregions 

0-  to  do  nothing 


IRA  energy  source  (for  impact  problems  ira  is  equal  to  0) 

ICM  1-  automatical  choice  of  time  step 

0 - time  step  =  abs(cmm)  *  stabl 


(cmm,stabl  see  below.) 


NBL  is  the  number  of  time  steps  to  do  before  copying  of  all 
flow  parameters  to  output  file. 

MI  1-  to  search  the  optimum  location  of  4  angle  tops  of 
rectangular  subregion  along  the  boundary 

ID1,ID2  are  the  numbers  of  time  steps  to  do  before  the 
data  record  for  restart(idl-  before  the  first  record) 

IK  is  the  maximum  number  of  time  steps  to  do. 

IGR  is  the  number  of  time  steps  to  do  before  printing  packed 
information  about  aU  subregion's  boundaries. 

ILEN  is  the  number  of  time  steps  to  do  before  the  next  graphic 
data  record. 

If  ilen  less  0  then  data  for  the  graphics  will  be 
written  in  time  step=abs(ilen*cmm)  in  nanoseconds. 

(cmm  parameter  described  below) 

PAR  is  the  number  of  flow  parameters  which  must  be  written  for 
graphics  data,  (from  5  to  7  parameters) 

rVSP  if  ivsp  is  not  equal  to  0  then  ivsp  is  input  parameter  for 
subroutine  AUTOSP,  which  changes  types  of "  contact " 
between  subregions  automatically  . 

This  automatical  regimes  must  be  programmed  by  user  for  the 
concrete  task  in  subroutine  AUTOSP. 

CMM  is  parameter  which  determines  time  step  (stt) 
stt  =  cmm  *  stabl  if  cmm>0. 

If  cmm<0.  then  cmm  will  be  read  from  restart  or 
cmm=  abs(cmm)  in  case  of  initial  data. 

H3  space  precision.(H3  equals  to  0.1  of  cell  size) 

H3  is  widely  used  in  the  code  as  a  characteristic  size. 

HDEL 

The  boundary  section  with  length  <  hdel*(cell  size)  will  be 
combined  with  neighboring  one. 

T  time  cpu  hmit  in  seconds 

STABL  stabiUty  coefficient  stt=  cmm  *  stabl 


The  second  part  of  the  initial  data  file  includes  an  information  about 
geometric  configurations  of  the  problem  to  be  solved.  The  information  may 


be  prepared  using  so  called  UNIVERSAL  INPUT.  In  this  case  the  first 
stringin  the  second  part  of  the  initial  data  file  must  be 
UNIVER . 

One  Universal  input  starting  with  UNIVER.. and  finishing  with 
UNIEND  specifies  only  one  subregion  (numerical  block). 

To  specify  several  subregions  one  must  use  this  block  of 
input  strings  several  times. 

The  boundaries  of  numerical  subregion  are  specified  in 
accordance  with  the  following  rules: 

1)  If  the  problem  contains  the  rigid  walls  and  shock  fronts 
they  must  be  specified  right  after  the  string  UNIVER.... 

For  example,  in  the  case  of  rigid  wall  the  input  looks  like 
this 

UNIVER . 

RIGIDRD 

rWALL,  INW,  ITP,  AW,  BW,  CW;STRAIGHT  LINE:F=INW*(AW*X+BW*Y-CW)=0. 

1  ,  -1  ,  1,  0.  ,  1.,  0. 

-1  ,  -1  ,  1,  0.  ,  1.,  0. 

Here  the  second  string  (after  UNIVER)  starting  as  RIGIDRD  indicates 
to  the  beginning  of  wall  specification,  which  is  finished  when  the  first 
parameter  (IWALL)  equals  to  some  negative  value  (-1  in  this  example).  In 
the  example  presented  above  the  rigid  wall  specified  is  the  straight  line 
(ITP=1)  and  input  parameters  are  the  coefficients  which  determine  it's 
orientation.  There  are  following  types  of  walls 
ITP=1 -rigid  wall  specified  by  the  straight  line 
ITP=2-rigid  wall  specified  by  the  circle 

ITP=3-  "transparent"  wall  (matter  can  flow  out  through  this  wall) 

To  specify  the  front  one  must  put  the  string  starting  as  FRNTREAD 
either  after  the  UNIVER...  or  after  the  rigid  wall  specification.  The  front 
specification  looks  hke  the  following 

FRNTREAD 

cl  —  read  of  front  information 
ifT,ityp,  1  the_number_of_front,  the_type_of_front 
icont,  I  the  number  of  contour 

a,b,c,  I  the  contour  is  straight  line  a*x+b*y-c=0  !  icont=nomco+ 1 
ibeg  end,!  l-begin,2-end 
-1 


Here  the  expressions  on  the  right  from  |  are  comments.  The  front 
input  is  finished  when  ifr  equals  to  some  negative  value  (-1  in  this 
example).  The  beginning  and  the  end  of  the  front  move  along  the  stright 
fines  (called  contour).  There  are  not  more  than  two  contours  corresponding 
to  the  same  front.  The  sequence  of  the  numbers  of  the  contours  specified 
must  be  increasing.  Every  contour  corresponds  either  to  the  beginning  of 


the  front  (ibeg_end=l)  or  to  the  end.  Using  such  an  input  several  fronts 
can  be  specified.  There  are  several  front  types  considered 
(the_type_of_front).  The  main  of  them  are 

1-  shock  front  in  undisturbed  soUd  media 

2- shock  front  in  undisturbed  porous  media 

3-  front  of  onedimensionahty 

Specification  of  aU  fronts  and  walls  must  be  closed  by  the  string 

END . 

2)  The  next  two  strings 

5.,  0.,  0.,::-USRT-,— VSTR— ,VPOL 
12,::-NUMBER  OF  URS 

give  the  velocity  components  (USRT,  VSTR),  porosity  of  the  subregion 
(VPOL)  and  number  of  material  (EOS).  Then  the  boundary  segments 
must  be  specified.  Boundary  orientation  is  counter-clockwise.  The  boundary 
segment  can  be  straight  hne,  circle  or  some  other  arbitrary  curve  specified 
in  the  table  form.  The  next  string  gives  the  coordinates  of  the  first  point 
of  the  boundary. 

0.,0.,::X0,Y0  -  the  location  of  the  first  point  for  a  boundary  (cm) 

If  the  first  boundary  segment  is  a  straight  hne,  then  the  next 
two  strings  are 

STRAIGHT  1 

0.,1.,  XE,YE,::  IS  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

If  the  next  segment  is  a  fragment  of  a  circle  then  must  be 
CIRCLE 

XCjYCjALF,::  is  the  center  of  the  CIRCLE(CM),ANGLE  (0-360  grad) 

The  boimdary  segment  can  be  any  arbitrary  curve.  In  this  case  it 
is  specified  using  the  table  input 

TABLE 

XE,YE,::  IS  THE  NODES  COORDINATES 

.......  ,  TABLE  INPUT  UNTIL  INPUT  ERROR  OCCURES 

END  OF  THE  TABLE 

The  specification  of  the  boundary  of  a  subregion  is  closed  by 
the  string 


REGEND 


The  four  angles  of  subregion  must  be  determined  for  the  purpose  of  grid 
generation,  when  the  boundary  of  the  subregion  is  specified.  The  angles 
input  looks  hke  this 

3— 4  ANGLES  LOCATION 

11,12,13,14,"  TO  FIND  4  ANGLES  OF  THE  RECTANGULAR  REGION  USING 
CORRESPONDENT  BOUNDARY  SEGMENTS  NUMBERS  11,12,13,14 

In  the  example  of  input  shown  below 

the  ANGLES  with  numbers  1,2,4  are  determined  automatically  and  the 
angle  3  is  found  in  the  beginning  of  the  boundary  segment  number  5. 
Initially  all  the  boimdary  segments  have  the  type  of  "FREE  BOUNDARY". 
If  one  needs  to  assign  some  other  boundary  type  to  the  boimdary  segment 
the  following  strings  are  necessary: 

BOUNDARY  TYPE  INPUT 
Then  depending  on  the  boundary  type  can  be: 

a) in  the  case  of  rigid  wall 

RIGID  ,;CHANGES  "FREE"  TYPE  BY  "RIGID" 

IREG,IRG  ,;  ireg-number  of  boundary  segment 

IRG-NUMBER  OF  RIGID  WALL 

b) in  the  case  of  shock  front 

FRONT  ,;CHANGES  "FREE"  TYPE  BY  "FRONT’ 

IREG,IFR  ,;  ireg-number  of  boundary  segment 

IFR-NUMBER  OF  FRONT 

c)  in  the  case  of  contact  discontinuity 
CONTACT  ,;  CHANGES  "free"  type  by  "contact’ 

IREG,N1,IREG1,;  IREG-NUMBER  OF  BOUNDARY  SEGMENT  FOR  CURRENT  SUBREGION,  Nl- 
NUMBER  OF  SUBREGION  WHICH  IS  IN  CONTACT,  IREGl-NUMBER  OF  BOUNDARY 
SEGMENT  FOR  N1  SUBREGION 

The  input  of  data  for  the  current  subregion  is  closed  by 
UNIEND 
Remark: 

We  have  described  UNIVERSAL  INPUT  in  the  text  above.  If  there  are 
special  program  for  more  simple  input  for  concrete  configuration  then  the 
name  of  this  program  is  the  first  string  in  the  second  part  of  the  input  data. 
For  example: 

DATSTl  ;  -name  of  the  program  which  treats  input  data 


How  to  run  the  code 


1)  Run  the  module  STARTT  in  the  directory  containing  file 
"NAM".DAT 

2)  Type  "NAM"  when  the  question  "ENTER  DATA  NAME"  appears 


module  STARTT  will  create  the  files: 

"NAM".TMP-  for  following  calculations 
"NAM"000.DAT-  first  file  for  graphics 
"NAM".CFG-  file  containing  some  parameters  (if  this  file  does  not  exist) 

An  example  of  CFG  file  is  presented  below 


- -VRB2.CFG - 

1.  ,  PLIM  [bar]  fine  1 

2.7104001E-05,  RUM  [G/(CM**3)]  Une  2 
13.93920  ,  EPSGAS  [M/C]  line3 

278.7840  ,  EPSCON  [M/C]  line4 

9.9999998E-03,  ALMIN  lineS 

0.3000000  ,  SINCOND  line6 

2.0000000E-02,  ALMN  hne7 

9.9999998E-03,  SHISAMI  Hne8 

5.000000  lE-02,  SHISAMA  line9 

0.7000000  ,  TET  HnelO 

0.000000  ,  ADP  hnell 


In  this  file  the  first  two  strings  are  minimal  pressure  and  density,  the  next 
two  strings  are  accuracy  of  calculation  of  wave  velocities  for  the  cases  of  gas 
EOS  and  table  EOS  correspondingly.  Parameters  ALMIN,  SINCOND 
(SINCOND  is  changing  from  0  to  1),  ALMN  are  responsible  for  the 
movement  of  boundary  nodes  on  the  boundary  of  "special"  type.  SHISAMI 
and  SHISAMA  are  minimum  and  maximum  shift  of  the  boundary  nodes 
during  the  procedure  of  boundary  node  distribution  (in  relative  units).  TET 
is  a  parameter  determining  the  convergence  of  grid  generation  (TET>0) 
and  ADP  is  a  level  of  grid  adaptation  to  solution.  In  the  present  version  of 
the  code  adaptation  in  not  included  (ADP=0) 

When  the  STARTT  finishes  its  work  the  following  message  appears 
THE  NUMBER  OF  PICTURE  =0 

3)  Run  the  module  BIG2T  (or  BIGIT)  to  continue  calculations 

4)  Reply  the  questions  appearing  on  the  screen 

To  avoid  the  grid  tangle  in  the  cases  of  large  deformations  to  split  the 
numerical  regions  into  subregions  is  recommended.  The  most  frequently 
appearing  question  concerning  that  is 


IF  YOU  DON'T  WANT  TO  SPLIT  SUBREGION  * 


ENTER  number  of  steps  without  splitting  analysis 

Enter  0  on  this  question  if  you  want  to  spht  subregion,  which  the  program 
suggests  to  spht.  Then  press  enter  to  continue  computations. 

The  question  ENTER  N,M  ACCORDING  TO  THE  STENCIL  means 
that  you  can  change  the  type  of  contact  between  subregions  N  and  M. 

Enter  -1,1  to  ignor  this  (and  -1,1  for  the  next  question) 

The  program  periodically  creates  the  files  "NAM"xxx.DAT  for  graphics 
and  renews  file  "NAM".TMP  for  restart. 

When  the  program  has  done  the  calculations  it  starts  the  dialogue 
from  the  beginning  (asking  "ENTER  DATA  NAME").  To  continue 
the  computations  of  the  same  problem  enter  the  same  name  NAM. 

5)  Analyze  the  results  of  computation  using  graphics 
modules.  For  this  purpose  created  graph  files 
"NAM"***.DAT  (and  k"NAM"***.DAT)  must  be  copied  to  IBM  PC 
Then  they  can  be  depicted  on  the  screen  using  and 

written  in  TIF  file  using  GR.EXE 
(See  the  manual  for  GR.EXE) 

6)  Continue  calculations  (point  3)  if  it  is  necessary 

7)  Change  configuration  of  the  task,if  it  is  necessary 
if  you  want  to  change  configuration  you  must  create 
the  file  D"NAM".DAT  (see  description  of  the  file) 
and  run  the  program  REP 

Brief  manual  for  GR.EXE 

Executing  module  GR.EXE  shows  the  results  of  computation  on  the 
screen  of  IBM  PC  (DOS  is  required)  and  writes  the  pictures  shown  on  the 
screen  in  TIF  format.  Besides  the  resulting  files  "NAME"***.DAT  to  be 
shown,  one  more  file  K"NAME".DAT is  necessary  to  determine  the  limits 
of  the  picture.  This  file  can  be  created  anytime  by  coping  already  existing 
files  of  such  type.  After  running  GR.EXE  one  need  to  answer  the  questions 
in  the  dialog. 

1)  Type  the  root  name  of  *.DAT  files  to  be  plotted  ("NAME") 

2)  Enter  the  number  of  picture  N  (file  "NAME'N.DAT  will  be  plotted) 

3)  Choose  the  flow  parameter  to  be  plotted  (1-density,  2-pressure,3- 
porosity,  4, 5-velocities,  6-energy) 

4)  Choose  the  type  of  picture  (I-levels  of  constant  values  and  numerical 
grid,  12-levels  of  constant  values  for  two  different  flow  parameters,  DY- 
cross  sections  along  the  axis  of  symmetry,  DX-  perpendicular  to  the  axis) 


5)  When  the  picture  is  shown  on  the  screen  one  can  plot  over  this  picture 
another  pictures  (if  the  axes  are  not  plotted  yet).  For  this  purpose  enter 
negative  number  -N  (where  N  is  the  number  of  picture  to  be  plotted) 

6)  To  Gnish  plotting  enter  any  character  (for  example  )  and  axes  will  be 
shown.  Then  type  enter ,  the  name  of  the  TIP  Gle  where  this  picture  will 
be  saved  and  CTRL-C  to  leave  the  program. 

7)  REMARK:  When  using  DY  or  DX  options  (plotting  cross  sections)  use 
CTRL-Z  to  return  back  for  plotting  pictures  aGer  the  question 

1=  INPUT  Y/X  COORDINATE  FOR  SECTION. 

There  is  an  example  of  GR.EXE  execution  presented  on  the  picture  below 
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Figure  .1  An  Example  of  GR.EXE  plot 

Example  of  initial  data  Gles 

Let  VIS  consider  some  demonstrative  examples  of  initial  data  generation 
for  hypervelocity  impact  problems 


EXAMPLE  1:  initial  data  file  VRB3.DAT 


Impact  of  W  projectile  on  A1  target 
JSPR-LOS-LST-IRD-ISP-IRA-ICM-AIR 
1  1  1  2  1  00  1  0 
NBL-MI-1D1-ID2-1K-1GR-1LENPARIVSP 
9900  1  30  30  300  30-500  6  0 

- CMM— : - H3 - : - HDEL— : - T - STABL— : 

-l.OOOOOE-03  0.02000E-00  0.06000E+00  l.lOOOOE+05  0.60000E+00 

UNIVER  -NAME  OF  PROGRAM  WHICH  TREATS  INPUT  DATA 
RIGIDRD  ,  ;  READING  AXIS  AND  RIGID  WALLS 

C  IWALL,  INW,  ITP,  AW ,  BW,  CW  ;STRAIGHT  LINE;  F=AW*X+BW*Y-CW=0. 

1  ,  -1  ,  1,  0.  ,  1.,  0. 


-1  ,  -1  ,  1,  0.  ,  1.,  0.  exit, A, A, A, A 

END 

-5.,  0.,  0.,::-USRT-,— VSTR— ,VPOL 

44,;;-.NUMBER  of  EOS 

0.,0.,::  xo.YO  -  the  location  of  the  first  point  for  the  boundary  (cm) 

STRAIGHT  1 

1.. 0.,  XE,YE,:;  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  2 

1.. 4.,  XE,YE,::  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  3 

0.,4.,  XE,YE,:;  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  4 

0.,0.,  XE,YE,::  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

REGEND 
c-l -2-3-4 
1,2, 3, 4, 

BOUNDARY  TYPE  INPUT 

RIGID  ,;CHANGING  "FREE"  TYPE  BY  "RIGID" 

1.1, ;  IREG-NUMBER  OF  BOUNDARY  SEGMENT, IRG-NUMBER  OF  RIGID  WALL 

UNIEND  1 

UNFVER  ,  -NAME  OF  PROGRAM  WHICH  TREATS  INPUT  DATA 

END 

0.,  0.,  0.,;:-USRT-,— VSTR— ,VPOL 

41,  ::--NUMBEROFEOS 

-8.5,0.,::  xo.yo  -  the  location  of  the  first  point  for  the  boundary  (cm) 

STRAIGHT  1 

0.,0.,  XE,YE,:;  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

TABLE  1 

0.,0., 

0.,8.5, 

0.,8.5, 

-8.5, 8.5, 

-8. 5, 8. 5, 

-8.5,0., 

-8.5,0., 

END  OF  TABLE 
REGEND 

3— 4 

1)2,„ 

BOUNDARY  TYPE  INPUT 
RIGID 

1.1, 

UNIEND  2 
ENDINPUT 


Fig.l  Initial  configuration  for  initial  data  file  VRB3.DAT 


EXAMPLE  2;  initial  data  file  VRB2.DAT 


IMPACT  with  shock  fitting 

JSPR-LOS-LST-IRD-ISP-IRA-ICM-AIR 
1  1  1  2  1  00  1  0 
NBL-MI-ID1-ID2-IK-IGR-ILENPARIVSP 

9900  1  30  30  300  30-500  6  0 

- CMM— : - H3 - ; - HDEL— : - T - STABL— : 

-l.OOOOOE-03  0.02000E-00  0.06000E+00  l.lOOOOE+05  0.60000E+00 

UNIVER  , -NAME  OF  PROGRAM  WHICH  TREATS  INPUT  DATA 

FRNTREAD 

Cl  —  read  of  the  front  information 

1,1,  ;c2  ifi:,ityp  ,:  the_number_of_front,  the_type_of_front 

1,  ;c3  icont,:  the_number_of_contour 

1.. 0..0.,  ;c4  a,b,c,:  the  contour  is  straight  line  a*x+b*y-c=0  ! 

icont=nomco+l 

1,  ;c5  ibeg_end,:  l-begin,2-end 

-1,  ;c6  end 

RIGIDRD  ,  5  READ  OF  AXIS  AND  RIGID  WALLS 

C  IWALL,  INW,  ITP,  AW  ,  BW,  CW  ;STRAIGHT  UNE:  F=AW*X+BW*Y-CW=0. 

1  ,  -1  ,  1,  0.  ,  1.,  0. 

-1  ,  -1  ,  1,  0.  ,  1.,  0.  exit,A,A,A,A 
END 

-5.,  0.,  0.,::-USRT-,— VSTR— ,VPOL 

44,  ::-NUMBER  OF  EOS 

0.,  0.,  xo.Yo  -  the  location  of  the  first  point  for  the  boundary  (cm) 
STRAIGHT  1 

1.. 0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  2 


l.,4.,  XE,YE,:;  ARE  COORDINATES  OF  THE  END  OF  THE 
STRAIGHT  SEGMENT  (CM) 

STRAIGHT  3 

0.,4.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  4 

0.,0.,  XE,YE,;:  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

REGEND 

1.. .2—3—4  ANGLE  LOCATION 
1, 2, 3,4, 

BOUNDARY  TYPE  INPUT 

RIGID  ,;YOU  WILL  CHANGE  "FREE"  TYPE  ON  "RIGID" 

1.1, ;  IREG-NUMBER  OF  BOUNDARY  SEGMENT, IRG-NUMBER  OF  RIGID  WALL 

UNIEND  1 

UNIVER  ,;  -NAME  OF  PROGRAM  WHICH  TREATS  INPUT  DATA 

END 

0.,  0.,  0.,::— USRT-,— VSTR— ,VPOL 

41,  "-NUMBER  OF  EOS 

-0.25,  0.,  X0,Y0  -  the  location  of  the  first  point  for  the  boundary  (CM) 

STRAIGHT  1 

0,,0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

STRAIGHT  2 

0.,4. 1,  XE,YE,;:  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

TABLE  3 

0.,4.1,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 
-0.,4.1,  XE,YE,;;  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

-0.25,3.0,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

-0.25,3.0,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

-0.25,0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

-0.25,0.,  XE,YE,::  ARE  COORDINATES  OF  THE  END  OF  THE  STRAIGHT  SEGMENT  (CM) 

END  OF  TABLE 
REGEND 

1.. .2,3—4  ANGLE  LOCATION 
lj2,3„ 

BOUNDARY  TYPE  INPUT 

RIGID  ,;YOU  WILL  CHANGE  "FREE"  TYPE  ON  "RIGID" 

1.1, ;  IREG-NUMBER  OF  BOUNDARY  SEGMENT,IRG-NUMBER  OF  RIGID  WALL 

FRONT 

3.1, ;  IREG-NUMBER  OF  BOUNDARY  SEGMENT, IRG-NUMBER  OF  FRONT 

UNIEND  2 
ENDINPUT 


Note:  Boundary  specified  by  TABLE  is  considered  as  one  boundary 
segment !!! 
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Fig.2  Initial  configuration  for  initial  data  file  VRB2.DAT 


Strategy  of  computations  and  some  useful  details 


Grid  splitting 

To  prevent  the  code  failure  in  grid  generation,  please,  avoid  numerical 
subregion  being  of  complicated  shape.  In  the  most  of  cases  when  it  is 
necessary,  the  program  analyses  the  form  of  regions  and  offers  the  ways  of 
dividing  them  into  subregions.  Nevertheless,  one  can  split  any  numerical 
region  running  module  REP.OUT.  This  module  requires  the  data  file 
D"name".DAT,  where  action  is  specified.  In  the  case  of  sphttmg  the  Erst 
two  strings  in  D"name".DAT must  be 
SPLSUB 
ENDINPUT 

Then  running  REP.  OUT  one  should  answer  the  following  equations: 

1)  Enter  the  number  of  the  region  to  be  split 

2)  Enter  option  to  reSne  the  grid  in  the  subregion  (0-no  reEnement,  1- 
reEnement  along  i-  index,  2-alongj  index,  3-  along  both) 

3)  Enter  option  to  make  the  gnd  rough  (answers  are  similar  as  for  2)) 

4  )  Enter  1-to  split  along  i  gnd  Ene,  2-to  split  along  j  gnd  Une,  0- 

no  splitting. 

If  you  put  1  or  2  on  the  question  4)  there  will  be  one  more  question: 

5)  Enter  the  number  of  i(or  j)  Une  along  which  the  region  wiU  be  spUt. 

To  leave  dialog  enter  any  negative  number  of  the  region. 


Removal  of  numerical  region  N 


Sometimes  it  is  desirable  to  delete  the  numerical  region  which  is  not 
interesting.  To  do  this  put  the  following  four  strings  into  the  file 
D"name".DAT: 

DELSUB 

N 

-1 

ENDINPUT 

where  N-  is  the  number  of  subregion  to  be  deleted. 

Then  run  REP.OUT again. 


Addition  of  a  new  numerical  subregion 

The  code  allows  one  to  include  into  consideration  additional  numerical 

regions  as  soon  as  it  is  required  (for  example  in  the  case  of  an  impact  on 

multilayered  spaced  target).  To  do  this  put  all  information  about  these 

regions  into  D"nam".dat  according  to  initial  data  format  (starting  from 

UNIVER 

and  ending  with 

ENDINPUT 

Then  run  REP.  OUT  and  a  new  TMP  file  will  be  generated.  After  that 
proceed  with  doing  computations  with  BIG2T.OUT. 


Changing  boundary  type. 

Sometimes  it  is  necessary  to  change  the  type  of  boundary  between 
subregions.  For  example,  newly  appearing  subregions  (after  splitting 
procedure)  have  so  called  "special"  type  of  boundaries  between  each  other. 
To  change  this  boundary  into  contact  type  is  possible  running  BIG2T.OUT 
.For  this  purpose: 

1)  Put  IRD=1  and  ISP=1  in  "NAMF'.DAT 

2)  Run  BIG2T.OUT and  enter  the  following  answers  on  the  questions: 
N1,N2 

N1,N2,M 

N2,N1,M 

-1,1 

-1,1 

Here  N1  and  N2  are  the  numbers  of  subregions  where  the  boundary  type  is 
to  be  changed  from  current  type  to  contact  (if  M=0)  or  to  "special" 
boundary  (if  M=4). 


What  does  it  mean  and  What  to  do  if 


Output  print:  ???CHANGE  CM???  -  one  of  the  mesh  rows  (or  columns) 
becomes  too  thin  to  satisfy  the  Courant  stabiUty  condition.  In  this  case 
program  blocks  this  row  with  the  neighboring  one  during  the  flow 
parameter  recalculation  to  the  next  time  step. 

Output  print:  DJCO(2)  -  This  output  usually  occurs  after  sphtting 
procedure.  One  should  stop  computations  and  repeat  it  once  again  by 
BIG2T  without  sphtting  of  subregion  at  the  same  moment.  One  can  do  it 
several  steps  later. 

Output  print:  CATCH  STOP  -  something  is  happened.  Press  enter  to 
continue  computations. 

GENERAL  RECOMMENDATIONS 

A.  Do  not  spht  subregions  without  necessity  !! 

B.  Changing  parameters  SHISAMI  and  SHISAMA  one  can  influence  the 
grid  generation.  You  can  do  the  grid  more  orthogonal  making  them  greater 
and  more  uniform  making  them  less.  If  you  have  some  problems  with  the 
grid,  for  example,  in  the  case  when  the  grid  has  strongly  different  cell 
dimensions  in  different  parts  of  the  numerical  region,  please,  decrease  these 
parameters  to  generate  more  uniform  grid. 

C.  Please,  look  pereodicaUy  at  your  results  using  graphics 

If  you  use  PC  computer  you  can  run  CINEMA.  EXE  in  the  directory 
where  are  your  graphic  files  NAME*. DAT  are. 


GOOD  LUCK 


